home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / meschach / meschach4.shar < prev    next >
Encoding:
Text File  |  1994-06-07  |  133.9 KB  |  5,327 lines  |  [TEXT/ttxt]
  1. # to unbundle, sh this file (in an empty directory)
  2. echo zmachine.c 1>&2
  3. sed >zmachine.c <<'//GO.SYSIN DD zmachine.c' 's/^-//'
  4. -
  5. -/**************************************************************************
  6. -**
  7. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  8. -**
  9. -**                 Meschach Library
  10. -** 
  11. -** This Meschach Library is provided "as is" without any express 
  12. -** or implied warranty of any kind with respect to this software. 
  13. -** In particular the authors shall not be liable for any direct, 
  14. -** indirect, special, incidental or consequential damages arising 
  15. -** in any way from use of the software.
  16. -** 
  17. -** Everyone is granted permission to copy, modify and redistribute this
  18. -** Meschach Library, provided:
  19. -**  1.  All copies contain this copyright notice.
  20. -**  2.  All modified copies shall carry a notice stating who
  21. -**      made the last modification and the date of such modification.
  22. -**  3.  No charge is made for this software or works derived from it.  
  23. -**      This clause shall not be construed as constraining other software
  24. -**      distributed on the same medium as this software, nor is a
  25. -**      distribution fee considered a charge.
  26. -**
  27. -***************************************************************************/
  28. -
  29. -
  30. -/*
  31. -  This file contains basic routines which are used by the functions
  32. -  involving complex vectors.
  33. -  These are the routines that should be modified in order to take
  34. -  full advantage of specialised architectures (pipelining, vector
  35. -  processors etc).
  36. -  */
  37. -static    char    *rcsid = "$Id: zmachine.c,v 1.1 1994/01/13 04:25:41 des Exp $";
  38. -
  39. -#include    <math.h>
  40. -#include    "machine.h"
  41. -#include        "zmatrix.h"
  42. -
  43. -
  44. -/* __zconj__ -- complex conjugate */
  45. -void    __zconj__(zp,len)
  46. -complex    *zp;
  47. -int    len;
  48. -{
  49. -    int        i;
  50. -
  51. -    for ( i = 0; i < len; i++ )
  52. -    zp[i].im = - zp[i].im;
  53. -}
  54. -
  55. -/* __zip__ -- inner product
  56. -    -- computes sum_i zp1[i].zp2[i] if flag == 0
  57. -            sum_i zp1[i]*.zp2[i] if flag != 0 */
  58. -complex    __zip__(zp1,zp2,len,flag)
  59. -complex    *zp1, *zp2;
  60. -int    flag, len;
  61. -{
  62. -    complex    sum;
  63. -    int        i;
  64. -
  65. -    sum.re = sum.im = 0.0;
  66. -    if ( flag )
  67. -    {
  68. -    for ( i = 0; i < len; i++ )
  69. -    {
  70. -        sum.re += zp1[i].re*zp2[i].re + zp1[i].im*zp2[i].im;
  71. -        sum.im += zp1[i].re*zp2[i].im - zp1[i].im*zp2[i].re;
  72. -    }
  73. -    }
  74. -    else
  75. -    {
  76. -    for ( i = 0; i < len; i++ )
  77. -    {
  78. -        sum.re += zp1[i].re*zp2[i].re - zp1[i].im*zp2[i].im;
  79. -        sum.im += zp1[i].re*zp2[i].im + zp1[i].im*zp2[i].re;
  80. -    }
  81. -    }
  82. -
  83. -    return sum;
  84. -}
  85. -
  86. -/* __zmltadd__ -- scalar multiply and add i.e. complex saxpy
  87. -    -- computes zp1[i] += s.zp2[i]  if flag == 0
  88. -    -- computes zp1[i] += s.zp2[i]* if flag != 0 */
  89. -void    __zmltadd__(zp1,zp2,s,len,flag)
  90. -complex    *zp1, *zp2, s;
  91. -int    flag, len;
  92. -{
  93. -    int        i;
  94. -    LongReal    t_re, t_im;
  95. -
  96. -    if ( ! flag )
  97. -    {
  98. -    for ( i = 0; i < len; i++ )
  99. -    {
  100. -        t_re = zp1[i].re + s.re*zp2[i].re - s.im*zp2[i].im;
  101. -        t_im = zp1[i].im + s.re*zp2[i].im + s.im*zp2[i].re;
  102. -        zp1[i].re = t_re;
  103. -        zp1[i].im = t_im;
  104. -    }
  105. -    }
  106. -    else
  107. -    {
  108. -    for ( i = 0; i < len; i++ )
  109. -    {
  110. -        t_re = zp1[i].re + s.re*zp2[i].re + s.im*zp2[i].im;
  111. -        t_im = zp1[i].im - s.re*zp2[i].im + s.im*zp2[i].re;
  112. -        zp1[i].re = t_re;
  113. -        zp1[i].im = t_im;
  114. -    }
  115. -    }
  116. -}
  117. -
  118. -/* __zmlt__ scalar complex multiply array c.f. sv_mlt() */
  119. -void    __zmlt__(zp,s,out,len)
  120. -complex    *zp, s, *out;
  121. -register int    len;
  122. -{
  123. -    int        i;
  124. -    LongReal    t_re, t_im;
  125. -
  126. -    for ( i = 0; i < len; i++ )
  127. -    {
  128. -    t_re = s.re*zp[i].re - s.im*zp[i].im;
  129. -    t_im = s.re*zp[i].im + s.im*zp[i].re;
  130. -    out[i].re = t_re;
  131. -    out[i].im = t_im;
  132. -    }
  133. -}
  134. -
  135. -/* __zadd__ -- add complex arrays c.f. v_add() */
  136. -void    __zadd__(zp1,zp2,out,len)
  137. -complex    *zp1, *zp2, *out;
  138. -int    len;
  139. -{
  140. -    int        i;
  141. -    for ( i = 0; i < len; i++ )
  142. -    {
  143. -    out[i].re = zp1[i].re + zp2[i].re;
  144. -    out[i].im = zp1[i].im + zp2[i].im;
  145. -    }
  146. -}
  147. -
  148. -/* __zsub__ -- subtract complex arrays c.f. v_sub() */
  149. -void    __zsub__(zp1,zp2,out,len)
  150. -complex    *zp1, *zp2, *out;
  151. -int    len;
  152. -{
  153. -    int        i;
  154. -    for ( i = 0; i < len; i++ )
  155. -    {
  156. -    out[i].re = zp1[i].re - zp2[i].re;
  157. -    out[i].im = zp1[i].im - zp2[i].im;
  158. -    }
  159. -}
  160. -
  161. -/* __zzero__ -- zeros an array of complex numbers */
  162. -void    __zzero__(zp,len)
  163. -complex    *zp;
  164. -int    len;
  165. -{
  166. -    /* if a Real precision zero is equivalent to a string of nulls */
  167. -    MEM_ZERO((char *)zp,len*sizeof(complex));
  168. -    /* else, need to zero the array entry by entry */
  169. -    /******************************
  170. -    while ( len-- )
  171. -    {
  172. -    zp->re = zp->im = 0.0;
  173. -    zp++;
  174. -    }
  175. -    ******************************/
  176. -}
  177. -
  178. //GO.SYSIN DD zmachine.c
  179. echo zcopy.c 1>&2
  180. sed >zcopy.c <<'//GO.SYSIN DD zcopy.c' 's/^-//'
  181. -
  182. -/**************************************************************************
  183. -**
  184. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  185. -**
  186. -**                 Meschach Library
  187. -** 
  188. -** This Meschach Library is provided "as is" without any express 
  189. -** or implied warranty of any kind with respect to this software. 
  190. -** In particular the authors shall not be liable for any direct, 
  191. -** indirect, special, incidental or consequential damages arising 
  192. -** in any way from use of the software.
  193. -** 
  194. -** Everyone is granted permission to copy, modify and redistribute this
  195. -** Meschach Library, provided:
  196. -**  1.  All copies contain this copyright notice.
  197. -**  2.  All modified copies shall carry a notice stating who
  198. -**      made the last modification and the date of such modification.
  199. -**  3.  No charge is made for this software or works derived from it.  
  200. -**      This clause shall not be construed as constraining other software
  201. -**      distributed on the same medium as this software, nor is a
  202. -**      distribution fee considered a charge.
  203. -**
  204. -***************************************************************************/
  205. -
  206. -
  207. -static    char    rcsid[] = "$Id: zcopy.c,v 1.1 1994/01/13 04:28:42 des Exp $";
  208. -#include    <stdio.h>
  209. -#include    "zmatrix.h"
  210. -
  211. -
  212. -
  213. -/* _zm_copy -- copies matrix into new area */
  214. -ZMAT    *_zm_copy(in,out,i0,j0)
  215. -ZMAT    *in,*out;
  216. -u_int    i0,j0;
  217. -{
  218. -    u_int    i /* ,j */;
  219. -
  220. -    if ( in==ZMNULL )
  221. -        error(E_NULL,"_zm_copy");
  222. -    if ( in==out )
  223. -        return (out);
  224. -    if ( out==ZMNULL || out->m < in->m || out->n < in->n )
  225. -        out = zm_resize(out,in->m,in->n);
  226. -
  227. -    for ( i=i0; i < in->m; i++ )
  228. -        MEM_COPY(&(in->me[i][j0]),&(out->me[i][j0]),
  229. -                (in->n - j0)*sizeof(complex));
  230. -        /* for ( j=j0; j < in->n; j++ )
  231. -            out->me[i][j] = in->me[i][j]; */
  232. -
  233. -    return (out);
  234. -}
  235. -
  236. -/* _zv_copy -- copies vector into new area */
  237. -ZVEC    *_zv_copy(in,out,i0)
  238. -ZVEC    *in,*out;
  239. -u_int    i0;
  240. -{
  241. -    /* u_int    i,j; */
  242. -
  243. -    if ( in==ZVNULL )
  244. -        error(E_NULL,"_zv_copy");
  245. -    if ( in==out )
  246. -        return (out);
  247. -    if ( out==ZVNULL || out->dim < in->dim )
  248. -        out = zv_resize(out,in->dim);
  249. -
  250. -    MEM_COPY(&(in->ve[i0]),&(out->ve[i0]),(in->dim - i0)*sizeof(complex));
  251. -    /* for ( i=i0; i < in->dim; i++ )
  252. -        out->ve[i] = in->ve[i]; */
  253. -
  254. -    return (out);
  255. -}
  256. -
  257. -
  258. -/*
  259. -    The z._move() routines are for moving blocks of memory around
  260. -    within Meschach data structures and for re-arranging matrices,
  261. -    vectors etc.
  262. -*/
  263. -
  264. -/* zm_move -- copies selected pieces of a matrix
  265. -    -- moves the m0 x n0 submatrix with top-left cor-ordinates (i0,j0)
  266. -       to the corresponding submatrix of out with top-left co-ordinates
  267. -       (i1,j1)
  268. -    -- out is resized (& created) if necessary */
  269. -ZMAT    *zm_move(in,i0,j0,m0,n0,out,i1,j1)
  270. -ZMAT    *in, *out;
  271. -int    i0, j0, m0, n0, i1, j1;
  272. -{
  273. -    int        i;
  274. -
  275. -    if ( ! in )
  276. -    error(E_NULL,"zm_move");
  277. -    if ( i0 < 0 || j0 < 0 || i1 < 0 || j1 < 0 || m0 < 0 || n0 < 0 ||
  278. -     i0+m0 > in->m || j0+n0 > in->n )
  279. -    error(E_BOUNDS,"zm_move");
  280. -
  281. -    if ( ! out )
  282. -    out = zm_resize(out,i1+m0,j1+n0);
  283. -    else if ( i1+m0 > out->m || j1+n0 > out->n )
  284. -    out = zm_resize(out,max(out->m,i1+m0),max(out->n,j1+n0));
  285. -
  286. -    for ( i = 0; i < m0; i++ )
  287. -    MEM_COPY(&(in->me[i0+i][j0]),&(out->me[i1+i][j1]),
  288. -         n0*sizeof(complex));
  289. -
  290. -    return out;
  291. -}
  292. -
  293. -/* zv_move -- copies selected pieces of a vector
  294. -    -- moves the length dim0 subvector with initial index i0
  295. -       to the corresponding subvector of out with initial index i1
  296. -    -- out is resized if necessary */
  297. -ZVEC    *zv_move(in,i0,dim0,out,i1)
  298. -ZVEC    *in, *out;
  299. -int    i0, dim0, i1;
  300. -{
  301. -    if ( ! in )
  302. -    error(E_NULL,"zv_move");
  303. -    if ( i0 < 0 || dim0 < 0 || i1 < 0 ||
  304. -     i0+dim0 > in->dim )
  305. -    error(E_BOUNDS,"zv_move");
  306. -
  307. -    if ( (! out) || i1+dim0 > out->dim )
  308. -    out = zv_resize(out,i1+dim0);
  309. -
  310. -    MEM_COPY(&(in->ve[i0]),&(out->ve[i1]),dim0*sizeof(complex));
  311. -
  312. -    return out;
  313. -}
  314. -
  315. -
  316. -/* zmv_move -- copies selected piece of matrix to a vector
  317. -    -- moves the m0 x n0 submatrix with top-left co-ordinate (i0,j0) to
  318. -       the subvector with initial index i1 (and length m0*n0)
  319. -    -- rows are copied contiguously
  320. -    -- out is resized if necessary */
  321. -ZVEC    *zmv_move(in,i0,j0,m0,n0,out,i1)
  322. -ZMAT    *in;
  323. -ZVEC    *out;
  324. -int    i0, j0, m0, n0, i1;
  325. -{
  326. -    int        dim1, i;
  327. -
  328. -    if ( ! in )
  329. -    error(E_NULL,"zmv_move");
  330. -    if ( i0 < 0 || j0 < 0 || m0 < 0 || n0 < 0 || i1 < 0 ||
  331. -     i0+m0 > in->m || j0+n0 > in->n )
  332. -    error(E_BOUNDS,"zmv_move");
  333. -
  334. -    dim1 = m0*n0;
  335. -    if ( (! out) || i1+dim1 > out->dim )
  336. -    out = zv_resize(out,i1+dim1);
  337. -
  338. -    for ( i = 0; i < m0; i++ )
  339. -    MEM_COPY(&(in->me[i0+i][j0]),&(out->ve[i1+i*n0]),n0*sizeof(complex));
  340. -
  341. -    return out;
  342. -}
  343. -
  344. -/* zvm_move -- copies selected piece of vector to a matrix
  345. -    -- moves the subvector with initial index i0 and length m1*n1 to
  346. -       the m1 x n1 submatrix with top-left co-ordinate (i1,j1)
  347. -        -- copying is done by rows
  348. -    -- out is resized if necessary */
  349. -ZMAT    *zvm_move(in,i0,out,i1,j1,m1,n1)
  350. -ZVEC    *in;
  351. -ZMAT    *out;
  352. -int    i0, i1, j1, m1, n1;
  353. -{
  354. -    int        dim0, i;
  355. -
  356. -    if ( ! in )
  357. -    error(E_NULL,"zvm_move");
  358. -    if ( i0 < 0 || i1 < 0 || j1 < 0 || m1 < 0 || n1 < 0 ||
  359. -     i0+m1*n1 > in->dim )
  360. -    error(E_BOUNDS,"zvm_move");
  361. -
  362. -    if ( ! out )
  363. -    out = zm_resize(out,i1+m1,j1+n1);
  364. -    else
  365. -    out = zm_resize(out,max(i1+m1,out->m),max(j1+n1,out->n));
  366. -
  367. -    dim0 = m1*n1;
  368. -    for ( i = 0; i < m1; i++ )
  369. -    MEM_COPY(&(in->ve[i0+i*n1]),&(out->me[i1+i][j1]),n1*sizeof(complex));
  370. -
  371. -    return out;
  372. -}
  373. //GO.SYSIN DD zcopy.c
  374. echo zmatio.c 1>&2
  375. sed >zmatio.c <<'//GO.SYSIN DD zmatio.c' 's/^-//'
  376. -
  377. -/**************************************************************************
  378. -**
  379. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  380. -**
  381. -**                 Meschach Library
  382. -** 
  383. -** This Meschach Library is provided "as is" without any express 
  384. -** or implied warranty of any kind with respect to this software. 
  385. -** In particular the authors shall not be liable for any direct, 
  386. -** indirect, special, incidental or consequential damages arising 
  387. -** in any way from use of the software.
  388. -** 
  389. -** Everyone is granted permission to copy, modify and redistribute this
  390. -** Meschach Library, provided:
  391. -**  1.  All copies contain this copyright notice.
  392. -**  2.  All modified copies shall carry a notice stating who
  393. -**      made the last modification and the date of such modification.
  394. -**  3.  No charge is made for this software or works derived from it.  
  395. -**      This clause shall not be construed as constraining other software
  396. -**      distributed on the same medium as this software, nor is a
  397. -**      distribution fee considered a charge.
  398. -**
  399. -***************************************************************************/
  400. -
  401. -
  402. -
  403. -#include        <stdio.h>
  404. -#include        <ctype.h>
  405. -#include        "zmatrix.h"
  406. -
  407. -static char rcsid[] = "$Id: zmatio.c,v 1.1 1994/01/13 04:25:18 des Exp $";
  408. -
  409. -
  410. -
  411. -/* local variables */
  412. -static char line[MAXLINE];
  413. -
  414. -/**************************************************************************
  415. -  Input routines
  416. -  **************************************************************************/
  417. -
  418. -complex    z_finput(fp)
  419. -FILE    *fp;
  420. -{
  421. -    int        io_code;
  422. -    complex    z;
  423. -
  424. -    skipjunk(fp);
  425. -    if ( isatty(fileno(fp)) )
  426. -    {
  427. -    do {
  428. -        fprintf(stderr,"real and imag parts: ");
  429. -        if ( fgets(line,MAXLINE,fp) == NULL )
  430. -        error(E_EOF,"z_finput");
  431. -#if REAL == DOUBLE
  432. -        io_code = sscanf(line,"%lf%lf",&z.re,&z.im);
  433. -#elif REAL == FLOAT
  434. -        io_code = sscanf(line,"%f%f",&z.re,&z.im);
  435. -#endif
  436. -
  437. -    } while ( io_code != 2 );
  438. -    }
  439. -    else
  440. -#if REAL == DOUBLE
  441. -      if ( (io_code=fscanf(fp," (%lf,%lf)",&z.re,&z.im)) < 2 )
  442. -#elif REAL == FLOAT
  443. -      if ( (io_code=fscanf(fp," (%f,%f)",&z.re,&z.im)) < 2 )
  444. -#endif
  445. -        error((io_code == EOF) ? E_EOF : E_FORMAT,"z_finput");
  446. -
  447. -    return z;
  448. -}
  449. -
  450. -
  451. -ZMAT    *zm_finput(fp,a)
  452. -FILE    *fp;
  453. -ZMAT    *a;
  454. -{
  455. -     ZMAT        *izm_finput(),*bzm_finput();
  456. -     
  457. -     if ( isatty(fileno(fp)) )
  458. -      return izm_finput(fp,a);
  459. -     else
  460. -      return bzm_finput(fp,a);
  461. -}
  462. -
  463. -/* izm_finput -- interactive input of matrix */
  464. -ZMAT     *izm_finput(fp,mat)
  465. -FILE    *fp;
  466. -ZMAT     *mat;
  467. -{
  468. -     char       c;
  469. -     u_int      i, j, m, n, dynamic;
  470. -     /* dynamic set to TRUE if memory allocated here */
  471. -     
  472. -     /* get matrix size */
  473. -     if ( mat != ZMNULL && mat->m<MAXDIM && mat->n<MAXDIM )
  474. -     {  m = mat->m;     n = mat->n;     dynamic = FALSE;        }
  475. -     else
  476. -     {
  477. -      dynamic = TRUE;
  478. -      do
  479. -      {
  480. -           fprintf(stderr,"ComplexMatrix: rows cols:");
  481. -           if ( fgets(line,MAXLINE,fp)==NULL )
  482. -            error(E_INPUT,"izm_finput");
  483. -      } while ( sscanf(line,"%u%u",&m,&n)<2 || m>MAXDIM || n>MAXDIM );
  484. -      mat = zm_get(m,n);
  485. -     }
  486. -     
  487. -     /* input elements */
  488. -     for ( i=0; i<m; i++ )
  489. -     {
  490. -     redo:
  491. -      fprintf(stderr,"row %u:\n",i);
  492. -      for ( j=0; j<n; j++ )
  493. -           do
  494. -           {
  495. -           redo2:
  496. -            fprintf(stderr,"entry (%u,%u): ",i,j);
  497. -            if ( !dynamic )
  498. -             fprintf(stderr,"old (%14.9g,%14.9g) new: ",
  499. -                 mat->me[i][j].re,mat->me[i][j].im);
  500. -            if ( fgets(line,MAXLINE,fp)==NULL )
  501. -             error(E_INPUT,"izm_finput");
  502. -            if ( (*line == 'b' || *line == 'B') && j > 0 )
  503. -            {   j--;    dynamic = FALSE;        goto redo2;     }
  504. -            if ( (*line == 'f' || *line == 'F') && j < n-1 )
  505. -            {   j++;    dynamic = FALSE;        goto redo2;     }
  506. -           } while ( *line=='\0' ||
  507. -#if REAL == DOUBLE
  508. -             sscanf(line,"%lf%lf",
  509. -#elif REAL == FLOAT
  510. -            sscanf(line,"%f%f",
  511. -#endif    
  512. -                &mat->me[i][j].re,&mat->me[i][j].im)<1 );
  513. -      fprintf(stderr,"Continue: ");
  514. -      fscanf(fp,"%c",&c);
  515. -      if ( c == 'n' || c == 'N' )
  516. -      {    dynamic = FALSE;                 goto redo;      }
  517. -      if ( (c == 'b' || c == 'B') /* && i > 0 */ )
  518. -      {     if ( i > 0 )
  519. -            i--;
  520. -        dynamic = FALSE;        goto redo;
  521. -      }
  522. -     }
  523. -     
  524. -     return (mat);
  525. -}
  526. -
  527. -/* bzm_finput -- batch-file input of matrix */
  528. -ZMAT     *bzm_finput(fp,mat)
  529. -FILE    *fp;
  530. -ZMAT     *mat;
  531. -{
  532. -     u_int      i,j,m,n,dummy;
  533. -     int        io_code;
  534. -     
  535. -     /* get dimension */
  536. -     skipjunk(fp);
  537. -     if ((io_code=fscanf(fp," ComplexMatrix: %u by %u",&m,&n)) < 2 ||
  538. -     m>MAXDIM || n>MAXDIM )
  539. -      error(io_code==EOF ? E_EOF : E_FORMAT,"bzm_finput");
  540. -     
  541. -     /* allocate memory if necessary */
  542. -     if ( mat==ZMNULL || mat->m<m || mat->n<n )
  543. -      mat = zm_resize(mat,m,n);
  544. -     
  545. -     /* get entries */
  546. -     for ( i=0; i<m; i++ )
  547. -     {
  548. -      skipjunk(fp);
  549. -      if ( fscanf(fp," row %u:",&dummy) < 1 )
  550. -           error(E_FORMAT,"bzm_finput");
  551. -      for ( j=0; j<n; j++ )
  552. -      {
  553. -          /* printf("bzm_finput: j = %d\n", j); */
  554. -#if REAL == DOUBLE
  555. -          if ((io_code=fscanf(fp," ( %lf , %lf )",
  556. -#elif REAL == FLOAT
  557. -          if ((io_code=fscanf(fp," ( %f , %f )",
  558. -#endif
  559. -                  &mat->me[i][j].re,&mat->me[i][j].im)) < 2 )
  560. -          error(io_code==EOF ? E_EOF : E_FORMAT,"bzm_finput");
  561. -      }
  562. -     }
  563. -     
  564. -     return (mat);
  565. -}
  566. -
  567. -ZVEC     *zv_finput(fp,x)
  568. -FILE    *fp;
  569. -ZVEC     *x;
  570. -{
  571. -     ZVEC        *izv_finput(),*bzv_finput();
  572. -     
  573. -     if ( isatty(fileno(fp)) )
  574. -      return izv_finput(fp,x);
  575. -     else
  576. -      return bzv_finput(fp,x);
  577. -}
  578. -
  579. -/* izv_finput -- interactive input of vector */
  580. -ZVEC     *izv_finput(fp,vec)
  581. -FILE    *fp;
  582. -ZVEC     *vec;
  583. -{
  584. -     u_int      i,dim,dynamic;  /* dynamic set if memory allocated here */
  585. -     
  586. -     /* get vector dimension */
  587. -     if ( vec != ZVNULL && vec->dim<MAXDIM )
  588. -     {  dim = vec->dim; dynamic = FALSE;        }
  589. -     else
  590. -     {
  591. -      dynamic = TRUE;
  592. -      do
  593. -      {
  594. -           fprintf(stderr,"ComplexVector: dim: ");
  595. -           if ( fgets(line,MAXLINE,fp)==NULL )
  596. -            error(E_INPUT,"izv_finput");
  597. -      } while ( sscanf(line,"%u",&dim)<1 || dim>MAXDIM );
  598. -      vec = zv_get(dim);
  599. -     }
  600. -     
  601. -     /* input elements */
  602. -     for ( i=0; i<dim; i++ )
  603. -      do
  604. -      {
  605. -      redo:
  606. -           fprintf(stderr,"entry %u: ",i);
  607. -           if ( !dynamic )
  608. -            fprintf(stderr,"old (%14.9g,%14.9g) new: ",
  609. -                vec->ve[i].re,vec->ve[i].im);
  610. -           if ( fgets(line,MAXLINE,fp)==NULL )
  611. -            error(E_INPUT,"izv_finput");
  612. -           if ( (*line == 'b' || *line == 'B') && i > 0 )
  613. -           {        i--;    dynamic = FALSE;        goto redo;         }
  614. -           if ( (*line == 'f' || *line == 'F') && i < dim-1 )
  615. -           {        i++;    dynamic = FALSE;        goto redo;         }
  616. -      } while ( *line=='\0' ||
  617. -#if REAL == DOUBLE
  618. -            sscanf(line,"%lf%lf",
  619. -#elif REAL == FLOAT
  620. -            sscanf(line,"%f%f",
  621. -#endif  
  622. -               &vec->ve[i].re,&vec->ve[i].im) < 2 );
  623. -     
  624. -     return (vec);
  625. -}
  626. -
  627. -/* bzv_finput -- batch-file input of vector */
  628. -ZVEC     *bzv_finput(fp,vec)
  629. -FILE    *fp;
  630. -ZVEC    *vec;
  631. -{
  632. -     u_int      i,dim;
  633. -     int        io_code;
  634. -     
  635. -     /* get dimension */
  636. -     skipjunk(fp);
  637. -     if ((io_code=fscanf(fp," ComplexVector: dim:%u",&dim)) < 1 ||
  638. -      dim>MAXDIM )
  639. -     error(io_code==EOF ? 7 : 6,"bzv_finput");
  640. -
  641. -     
  642. -     /* allocate memory if necessary */
  643. -     if ( vec==ZVNULL || vec->dim<dim )
  644. -      vec = zv_resize(vec,dim);
  645. -     
  646. -     /* get entries */
  647. -     skipjunk(fp);
  648. -     for ( i=0; i<dim; i++ )
  649. -#if REAL == DOUBLE
  650. -      if ((io_code=fscanf(fp," (%lf,%lf)",
  651. -#elif REAL == FLOAT
  652. -          if ((io_code=fscanf(fp," (%f,%f)",
  653. -#endif
  654. -                  &vec->ve[i].re,&vec->ve[i].im)) < 2 )
  655. -           error(io_code==EOF ? 7 : 6,"bzv_finput");
  656. -     
  657. -     return (vec);
  658. -}
  659. -
  660. -/**************************************************************************
  661. -  Output routines
  662. -  **************************************************************************/
  663. -static char    *zformat = " (%14.9g, %14.9g) ";
  664. -
  665. -char    *setzformat(f_string)
  666. -char    *f_string;
  667. -{
  668. -    char    *old_f_string;
  669. -    old_f_string = zformat;
  670. -    if ( f_string != (char *)NULL && *f_string != '\0' )
  671. -    zformat = f_string;
  672. -
  673. -    return old_f_string;
  674. -}
  675. -
  676. -void    z_foutput(fp,z)
  677. -FILE    *fp;
  678. -complex    z;
  679. -{
  680. -    fprintf(fp,zformat,z.re,z.im);
  681. -    putc('\n',fp);
  682. -}
  683. -
  684. -void    zm_foutput(fp,a)
  685. -FILE    *fp;
  686. -ZMAT     *a;
  687. -{
  688. -     u_int      i, j, tmp;
  689. -     
  690. -     if ( a == ZMNULL )
  691. -     {  fprintf(fp,"ComplexMatrix: NULL\n");   return;         }
  692. -     fprintf(fp,"ComplexMatrix: %d by %d\n",a->m,a->n);
  693. -     if ( a->me == (complex **)NULL )
  694. -     {  fprintf(fp,"NULL\n");           return;         }
  695. -     for ( i=0; i<a->m; i++ )   /* for each row... */
  696. -     {
  697. -      fprintf(fp,"row %u: ",i);
  698. -      for ( j=0, tmp=1; j<a->n; j++, tmp++ )
  699. -      {             /* for each col in row... */
  700. -           fprintf(fp,zformat,a->me[i][j].re,a->me[i][j].im);
  701. -           if ( ! (tmp % 2) )       putc('\n',fp);
  702. -      }
  703. -      if ( tmp % 2 != 1 )   putc('\n',fp);
  704. -     }
  705. -}
  706. -
  707. -void    zv_foutput(fp,x)
  708. -FILE    *fp;
  709. -ZVEC     *x;
  710. -{
  711. -     u_int      i, tmp;
  712. -     
  713. -     if ( x == ZVNULL )
  714. -     {  fprintf(fp,"ComplexVector: NULL\n");   return;         }
  715. -     fprintf(fp,"ComplexVector: dim: %d\n",x->dim);
  716. -     if ( x->ve == (complex *)NULL )
  717. -     {  fprintf(fp,"NULL\n");   return;         }
  718. -     for ( i=0, tmp=0; i<x->dim; i++, tmp++ )
  719. -     {
  720. -      fprintf(fp,zformat,x->ve[i].re,x->ve[i].im);
  721. -      if ( (tmp % 2) == 1 )   putc('\n',fp);
  722. -     }
  723. -     if ( (tmp % 2) != 0 )        putc('\n',fp);
  724. -}
  725. -
  726. -
  727. -void    zm_dump(fp,a)
  728. -FILE    *fp;
  729. -ZMAT     *a;
  730. -{
  731. -    u_int   i, j, tmp;
  732. -     
  733. -     if ( a == ZMNULL )
  734. -     {  fprintf(fp,"ComplexMatrix: NULL\n");   return;         }
  735. -     fprintf(fp,"ComplexMatrix: %d by %d @ 0x%lx\n",a->m,a->n,(long)a);
  736. -     fprintf(fp,"\tmax_m = %d, max_n = %d, max_size = %d\n",
  737. -         a->max_m, a->max_n, a->max_size);
  738. -     if ( a->me == (complex **)NULL )
  739. -     {  fprintf(fp,"NULL\n");           return;         }
  740. -     fprintf(fp,"a->me @ 0x%lx\n",(long)(a->me));
  741. -     fprintf(fp,"a->base @ 0x%lx\n",(long)(a->base));
  742. -     for ( i=0; i<a->m; i++ )   /* for each row... */
  743. -     {
  744. -      fprintf(fp,"row %u: @ 0x%lx ",i,(long)(a->me[i]));
  745. -      for ( j=0, tmp=1; j<a->n; j++, tmp++ )
  746. -      {             /* for each col in row... */
  747. -           fprintf(fp,zformat,a->me[i][j].re,a->me[i][j].im);
  748. -           if ( ! (tmp % 2) )       putc('\n',fp);
  749. -      }
  750. -      if ( tmp % 2 != 1 )   putc('\n',fp);
  751. -     }
  752. -}
  753. -
  754. -
  755. -
  756. -void    zv_dump(fp,x)
  757. -FILE    *fp;
  758. -ZVEC     *x;
  759. -{
  760. -     u_int      i, tmp;
  761. -     
  762. -     if ( ! x )
  763. -     {  fprintf(fp,"ComplexVector: NULL\n");   return;         }
  764. -     fprintf(fp,"ComplexVector: dim: %d @ 0x%lx\n",x->dim,(long)(x));
  765. -     if ( ! x->ve )
  766. -     {  fprintf(fp,"NULL\n");   return;         }
  767. -     fprintf(fp,"x->ve @ 0x%lx\n",(long)(x->ve));
  768. -     for ( i=0, tmp=0; i<x->dim; i++, tmp++ )
  769. -     {
  770. -      fprintf(fp,zformat,x->ve[i].re,x->ve[i].im);
  771. -      if ( tmp % 2 == 1 )   putc('\n',fp);
  772. -     }
  773. -     if ( tmp % 2 != 0 )        putc('\n',fp);
  774. -}
  775. -
  776. //GO.SYSIN DD zmatio.c
  777. echo zmemory.c 1>&2
  778. sed >zmemory.c <<'//GO.SYSIN DD zmemory.c' 's/^-//'
  779. -
  780. -/**************************************************************************
  781. -**
  782. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  783. -**
  784. -**                 Meschach Library
  785. -** 
  786. -** This Meschach Library is provided "as is" without any express 
  787. -** or implied warranty of any kind with respect to this software. 
  788. -** In particular the authors shall not be liable for any direct, 
  789. -** indirect, special, incidental or consequential damages arising 
  790. -** in any way from use of the software.
  791. -** 
  792. -** Everyone is granted permission to copy, modify and redistribute this
  793. -** Meschach Library, provided:
  794. -**  1.  All copies contain this copyright notice.
  795. -**  2.  All modified copies shall carry a notice stating who
  796. -**      made the last modification and the date of such modification.
  797. -**  3.  No charge is made for this software or works derived from it.  
  798. -**      This clause shall not be construed as constraining other software
  799. -**      distributed on the same medium as this software, nor is a
  800. -**      distribution fee considered a charge.
  801. -**
  802. -***************************************************************************/
  803. -
  804. -
  805. -/* Memory allocation and de-allocation for complex matrices and vectors */
  806. -
  807. -#include    <stdio.h>
  808. -#include    "zmatrix.h"
  809. -
  810. -static    char    rcsid[] = "$Id: zmemory.c,v 1.2 1994/04/05 02:13:14 des Exp $";
  811. -
  812. -
  813. -
  814. -/* zv_zero -- zeros all entries of a complex vector
  815. -   -- uses __zzero__() */
  816. -ZVEC    *zv_zero(x)
  817. -ZVEC    *x;
  818. -{
  819. -   if ( ! x )
  820. -     error(E_NULL,"zv_zero");
  821. -   __zzero__(x->ve,x->dim);
  822. -   
  823. -   return x;
  824. -}
  825. -
  826. -/* zm_zero -- zeros all entries of a complex matrix
  827. -   -- uses __zzero__() */
  828. -ZMAT    *zm_zero(A)
  829. -ZMAT    *A;
  830. -{
  831. -   int        i;
  832. -   
  833. -   if ( ! A )
  834. -     error(E_NULL,"zm_zero");
  835. -   for ( i = 0; i < A->m; i++ )
  836. -     __zzero__(A->me[i],A->n);
  837. -   
  838. -   return A;
  839. -}
  840. -
  841. -/* zm_get -- gets an mxn complex matrix (in ZMAT form) */
  842. -ZMAT    *zm_get(m,n)
  843. -int    m,n;
  844. -{
  845. -   ZMAT    *matrix;
  846. -   u_int    i;
  847. -   
  848. -   if (m < 0 || n < 0)
  849. -     error(E_NEG,"zm_get");
  850. -
  851. -   if ((matrix=NEW(ZMAT)) == (ZMAT *)NULL )
  852. -     error(E_MEM,"zm_get");
  853. -   else if (mem_info_is_on()) {
  854. -      mem_bytes(TYPE_ZMAT,0,sizeof(ZMAT));
  855. -      mem_numvar(TYPE_ZMAT,1);
  856. -   }
  857. -   
  858. -   matrix->m = m;        matrix->n = matrix->max_n = n;
  859. -   matrix->max_m = m;    matrix->max_size = m*n;
  860. -#ifndef SEGMENTED
  861. -   if ((matrix->base = NEW_A(m*n,complex)) == (complex *)NULL )
  862. -   {
  863. -      free(matrix);
  864. -      error(E_MEM,"zm_get");
  865. -   }
  866. -   else if (mem_info_is_on()) {
  867. -      mem_bytes(TYPE_ZMAT,0,m*n*sizeof(complex));
  868. -   }
  869. -#else
  870. -   matrix->base = (complex *)NULL;
  871. -#endif
  872. -   if ((matrix->me = (complex **)calloc(m,sizeof(complex *))) == 
  873. -       (complex **)NULL )
  874. -   {    free(matrix->base);    free(matrix);
  875. -    error(E_MEM,"zm_get");
  876. -     }
  877. -   else if (mem_info_is_on()) {
  878. -      mem_bytes(TYPE_ZMAT,0,m*sizeof(complex *));
  879. -   }
  880. -#ifndef SEGMENTED
  881. -   /* set up pointers */
  882. -   for ( i=0; i<m; i++ )
  883. -     matrix->me[i] = &(matrix->base[i*n]);
  884. -#else
  885. -   for ( i = 0; i < m; i++ )
  886. -     if ( (matrix->me[i]=NEW_A(n,complex)) == (complex *)NULL )
  887. -       error(E_MEM,"zm_get");
  888. -     else if (mem_info_is_on()) {
  889. -    mem_bytes(TYPE_ZMAT,0,n*sizeof(complex));
  890. -     }
  891. -#endif
  892. -   
  893. -   return (matrix);
  894. -}
  895. -
  896. -
  897. -/* zv_get -- gets a ZVEC of dimension 'dim'
  898. -   -- Note: initialized to zero */
  899. -ZVEC    *zv_get(size)
  900. -int    size;
  901. -{
  902. -   ZVEC    *vector;
  903. -
  904. -   if (size < 0)
  905. -     error(E_NEG,"zv_get");
  906. -
  907. -   if ((vector=NEW(ZVEC)) == (ZVEC *)NULL )
  908. -     error(E_MEM,"zv_get");
  909. -   else if (mem_info_is_on()) {
  910. -      mem_bytes(TYPE_ZVEC,0,sizeof(ZVEC));
  911. -      mem_numvar(TYPE_ZVEC,1);
  912. -   }
  913. -   vector->dim = vector->max_dim = size;
  914. -   if ((vector->ve=NEW_A(size,complex)) == (complex *)NULL )
  915. -   {
  916. -      free(vector);
  917. -      error(E_MEM,"zv_get");
  918. -   }
  919. -   else if (mem_info_is_on()) {
  920. -      mem_bytes(TYPE_ZVEC,0,size*sizeof(complex));
  921. -   }
  922. -   return (vector);
  923. -}
  924. -
  925. -/* zm_free -- returns ZMAT & asoociated memory back to memory heap */
  926. -int    zm_free(mat)
  927. -ZMAT    *mat;
  928. -{
  929. -#ifdef SEGMENTED
  930. -   int    i;
  931. -#endif
  932. -   
  933. -   if ( mat==(ZMAT *)NULL || (int)(mat->m) < 0 ||
  934. -       (int)(mat->n) < 0 )
  935. -     /* don't trust it */
  936. -     return (-1);
  937. -   
  938. -#ifndef SEGMENTED
  939. -   if ( mat->base != (complex *)NULL ) {
  940. -      if (mem_info_is_on()) {
  941. -     mem_bytes(TYPE_ZMAT,mat->max_m*mat->max_n*sizeof(complex),0);
  942. -      }       
  943. -      free((char *)(mat->base));
  944. -   }
  945. -#else
  946. -   for ( i = 0; i < mat->max_m; i++ )
  947. -     if ( mat->me[i] != (complex *)NULL ) {
  948. -    if (mem_info_is_on()) {
  949. -       mem_bytes(TYPE_ZMAT,mat->max_n*sizeof(complex),0);
  950. -    }
  951. -    free((char *)(mat->me[i]));
  952. -     }
  953. -#endif
  954. -   if ( mat->me != (complex **)NULL ) {
  955. -      if (mem_info_is_on()) {
  956. -     mem_bytes(TYPE_ZMAT,mat->max_m*sizeof(complex *),0);
  957. -      }       
  958. -      free((char *)(mat->me));
  959. -   }
  960. -   
  961. -   if (mem_info_is_on()) {
  962. -      mem_bytes(TYPE_ZMAT,sizeof(ZMAT),0);
  963. -      mem_numvar(TYPE_ZMAT,-1);
  964. -   }
  965. -   free((char *)mat);
  966. -   
  967. -   return (0);
  968. -}
  969. -
  970. -
  971. -/* zv_free -- returns ZVEC & asoociated memory back to memory heap */
  972. -int    zv_free(vec)
  973. -ZVEC    *vec;
  974. -{
  975. -   if ( vec==(ZVEC *)NULL || (int)(vec->dim) < 0 )
  976. -     /* don't trust it */
  977. -     return (-1);
  978. -   
  979. -   if ( vec->ve == (complex *)NULL ) {
  980. -      if (mem_info_is_on()) {
  981. -     mem_bytes(TYPE_ZVEC,sizeof(ZVEC),0);
  982. -     mem_numvar(TYPE_ZVEC,-1);
  983. -      }
  984. -      free((char *)vec);
  985. -   }
  986. -   else
  987. -   {
  988. -      if (mem_info_is_on()) {
  989. -     mem_bytes(TYPE_ZVEC,vec->max_dim*sizeof(complex)+
  990. -              sizeof(ZVEC),0);
  991. -     mem_numvar(TYPE_ZVEC,-1);
  992. -      }
  993. -      
  994. -      free((char *)vec->ve);
  995. -      free((char *)vec);
  996. -   }
  997. -   
  998. -   return (0);
  999. -}
  1000. -
  1001. -
  1002. -/* zm_resize -- returns the matrix A of size new_m x new_n; A is zeroed
  1003. -   -- if A == NULL on entry then the effect is equivalent to m_get() */
  1004. -ZMAT    *zm_resize(A,new_m,new_n)
  1005. -ZMAT    *A;
  1006. -int    new_m, new_n;
  1007. -{
  1008. -   u_int    i, new_max_m, new_max_n, new_size, old_m, old_n;
  1009. -   
  1010. -   if (new_m < 0 || new_n < 0)
  1011. -     error(E_NEG,"zm_resize");
  1012. -
  1013. -   if ( ! A )
  1014. -     return zm_get(new_m,new_n);
  1015. -   
  1016. -   if (new_m == A->m && new_n == A->n)
  1017. -     return A;
  1018. -
  1019. -   old_m = A->m;    old_n = A->n;
  1020. -   if ( new_m > A->max_m )
  1021. -   {    /* re-allocate A->me */
  1022. -      if (mem_info_is_on()) {
  1023. -     mem_bytes(TYPE_ZMAT,A->max_m*sizeof(complex *),
  1024. -              new_m*sizeof(complex *));
  1025. -      }
  1026. -
  1027. -      A->me = RENEW(A->me,new_m,complex *);
  1028. -      if ( ! A->me )
  1029. -    error(E_MEM,"zm_resize");
  1030. -   }
  1031. -   new_max_m = max(new_m,A->max_m);
  1032. -   new_max_n = max(new_n,A->max_n);
  1033. -   
  1034. -#ifndef SEGMENTED
  1035. -   new_size = new_max_m*new_max_n;
  1036. -   if ( new_size > A->max_size )
  1037. -   {    /* re-allocate A->base */
  1038. -      if (mem_info_is_on()) {
  1039. -     mem_bytes(TYPE_ZMAT,A->max_m*A->max_n*sizeof(complex),
  1040. -        new_size*sizeof(complex));      
  1041. -      }
  1042. -
  1043. -      A->base = RENEW(A->base,new_size,complex);
  1044. -      if ( ! A->base )
  1045. -    error(E_MEM,"zm_resize");
  1046. -      A->max_size = new_size;
  1047. -   }
  1048. -   
  1049. -   /* now set up A->me[i] */
  1050. -   for ( i = 0; i < new_m; i++ )
  1051. -     A->me[i] = &(A->base[i*new_n]);
  1052. -   
  1053. -   /* now shift data in matrix */
  1054. -   if ( old_n > new_n )
  1055. -   {
  1056. -      for ( i = 1; i < min(old_m,new_m); i++ )
  1057. -    MEM_COPY((char *)&(A->base[i*old_n]),
  1058. -         (char *)&(A->base[i*new_n]),
  1059. -         sizeof(complex)*new_n);
  1060. -   }
  1061. -   else if ( old_n < new_n )
  1062. -   {
  1063. -      for ( i = min(old_m,new_m)-1; i > 0; i-- )
  1064. -      {   /* copy & then zero extra space */
  1065. -     MEM_COPY((char *)&(A->base[i*old_n]),
  1066. -          (char *)&(A->base[i*new_n]),
  1067. -          sizeof(complex)*old_n);
  1068. -     __zzero__(&(A->base[i*new_n+old_n]),(new_n-old_n));
  1069. -      }
  1070. -      __zzero__(&(A->base[old_n]),(new_n-old_n));
  1071. -      A->max_n = new_n;
  1072. -   }
  1073. -   /* zero out the new rows.. */
  1074. -   for ( i = old_m; i < new_m; i++ )
  1075. -     __zzero__(&(A->base[i*new_n]),new_n);
  1076. -#else
  1077. -   if ( A->max_n < new_n )
  1078. -   {
  1079. -      complex    *tmp;
  1080. -      
  1081. -      for ( i = 0; i < A->max_m; i++ )
  1082. -      {
  1083. -     if (mem_info_is_on()) {
  1084. -        mem_bytes(TYPE_ZMAT,A->max_n*sizeof(complex),
  1085. -             new_max_n*sizeof(complex));
  1086. -     }
  1087. -
  1088. -     if ( (tmp = RENEW(A->me[i],new_max_n,complex)) == NULL )
  1089. -       error(E_MEM,"zm_resize");
  1090. -     else {
  1091. -        A->me[i] = tmp;
  1092. -     }
  1093. -      }
  1094. -      for ( i = A->max_m; i < new_max_m; i++ )
  1095. -      {
  1096. -     if ( (tmp = NEW_A(new_max_n,complex)) == NULL )
  1097. -       error(E_MEM,"zm_resize");
  1098. -     else {
  1099. -        A->me[i] = tmp;
  1100. -        if (mem_info_is_on()) {
  1101. -           mem_bytes(TYPE_ZMAT,0,new_max_n*sizeof(complex));
  1102. -        }
  1103. -     }
  1104. -      }
  1105. -   }
  1106. -   else if ( A->max_m < new_m )
  1107. -   {
  1108. -      for ( i = A->max_m; i < new_m; i++ )
  1109. -    if ( (A->me[i] = NEW_A(new_max_n,complex)) == NULL )
  1110. -      error(E_MEM,"zm_resize");
  1111. -    else if (mem_info_is_on()) {
  1112. -       mem_bytes(TYPE_ZMAT,0,new_max*sizeof(complex));
  1113. -    }
  1114. -      
  1115. -   }
  1116. -   
  1117. -   if ( old_n < new_n )
  1118. -   {
  1119. -      for ( i = 0; i < old_m; i++ )
  1120. -    __zzero__(&(A->me[i][old_n]),new_n-old_n);
  1121. -   }
  1122. -   
  1123. -   /* zero out the new rows.. */
  1124. -   for ( i = old_m; i < new_m; i++ )
  1125. -     __zzero__(A->me[i],new_n);
  1126. -#endif
  1127. -   
  1128. -   A->max_m = new_max_m;
  1129. -   A->max_n = new_max_n;
  1130. -   A->max_size = A->max_m*A->max_n;
  1131. -   A->m = new_m;    A->n = new_n;
  1132. -   
  1133. -   return A;
  1134. -}
  1135. -
  1136. -
  1137. -/* zv_resize -- returns the (complex) vector x with dim new_dim
  1138. -   -- x is set to the zero vector */
  1139. -ZVEC    *zv_resize(x,new_dim)
  1140. -ZVEC    *x;
  1141. -int    new_dim;
  1142. -{
  1143. -   if (new_dim < 0)
  1144. -     error(E_NEG,"zv_resize");
  1145. -
  1146. -   if ( ! x )
  1147. -     return zv_get(new_dim);
  1148. -
  1149. -   if (new_dim == x->dim)
  1150. -     return x;
  1151. -
  1152. -   if ( x->max_dim == 0 )    /* assume that it's from sub_zvec */
  1153. -     return zv_get(new_dim);
  1154. -   
  1155. -   if ( new_dim > x->max_dim )
  1156. -   {
  1157. -      if (mem_info_is_on()) { 
  1158. -     mem_bytes(TYPE_ZVEC,x->max_dim*sizeof(complex),
  1159. -              new_dim*sizeof(complex));
  1160. -      }
  1161. -
  1162. -      x->ve = RENEW(x->ve,new_dim,complex);
  1163. -      if ( ! x->ve )
  1164. -    error(E_MEM,"zv_resize");
  1165. -      x->max_dim = new_dim;
  1166. -   }
  1167. -   
  1168. -   if ( new_dim > x->dim )
  1169. -     __zzero__(&(x->ve[x->dim]),new_dim - x->dim);
  1170. -   x->dim = new_dim;
  1171. -   
  1172. -   return x;
  1173. -}
  1174. -
  1175. -
  1176. -/* varying arguments */
  1177. -
  1178. -#ifdef ANSI_C
  1179. -
  1180. -#include <stdarg.h>
  1181. -
  1182. -
  1183. -/* To allocate memory to many arguments. 
  1184. -   The function should be called:
  1185. -   zv_get_vars(dim,&x,&y,&z,...,NULL);
  1186. -   where 
  1187. -     int dim;
  1188. -     ZVEC *x, *y, *z,...;
  1189. -     The last argument should be NULL ! 
  1190. -     dim is the length of vectors x,y,z,...
  1191. -     returned value is equal to the number of allocated variables
  1192. -     Other gec_... functions are similar.
  1193. -*/
  1194. -
  1195. -int zv_get_vars(int dim,...) 
  1196. -{
  1197. -   va_list ap;
  1198. -   int i=0;
  1199. -   ZVEC **par;
  1200. -   
  1201. -   va_start(ap, dim);
  1202. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1203. -      *par = zv_get(dim);
  1204. -      i++;
  1205. -   } 
  1206. -
  1207. -   va_end(ap);
  1208. -   return i;
  1209. -}
  1210. -
  1211. -
  1212. -
  1213. -int zm_get_vars(int m,int n,...) 
  1214. -{
  1215. -   va_list ap;
  1216. -   int i=0;
  1217. -   ZMAT **par;
  1218. -   
  1219. -   va_start(ap, n);
  1220. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1221. -      *par = zm_get(m,n);
  1222. -      i++;
  1223. -   } 
  1224. -
  1225. -   va_end(ap);
  1226. -   return i;
  1227. -}
  1228. -
  1229. -
  1230. -
  1231. -/* To resize memory for many arguments. 
  1232. -   The function should be called:
  1233. -   v_resize_vars(new_dim,&x,&y,&z,...,NULL);
  1234. -   where 
  1235. -     int new_dim;
  1236. -     ZVEC *x, *y, *z,...;
  1237. -     The last argument should be NULL ! 
  1238. -     rdim is the resized length of vectors x,y,z,...
  1239. -     returned value is equal to the number of allocated variables.
  1240. -     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  1241. -     argument. 
  1242. -     Other *_resize_list() functions are similar.
  1243. -*/
  1244. -
  1245. -int zv_resize_vars(int new_dim,...)
  1246. -{
  1247. -   va_list ap;
  1248. -   int i=0;
  1249. -   ZVEC **par;
  1250. -   
  1251. -   va_start(ap, new_dim);
  1252. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1253. -      *par = zv_resize(*par,new_dim);
  1254. -      i++;
  1255. -   } 
  1256. -
  1257. -   va_end(ap);
  1258. -   return i;
  1259. -}
  1260. -
  1261. -
  1262. -
  1263. -int zm_resize_vars(int m,int n,...) 
  1264. -{
  1265. -   va_list ap;
  1266. -   int i=0;
  1267. -   ZMAT **par;
  1268. -   
  1269. -   va_start(ap, n);
  1270. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1271. -      *par = zm_resize(*par,m,n);
  1272. -      i++;
  1273. -   } 
  1274. -
  1275. -   va_end(ap);
  1276. -   return i;
  1277. -}
  1278. -
  1279. -
  1280. -/* To deallocate memory for many arguments. 
  1281. -   The function should be called:
  1282. -   v_free_vars(&x,&y,&z,...,NULL);
  1283. -   where 
  1284. -     ZVEC *x, *y, *z,...;
  1285. -     The last argument should be NULL ! 
  1286. -     There must be at least one not NULL argument.
  1287. -     returned value is equal to the number of allocated variables.
  1288. -     Returned value of x,y,z,.. is VNULL.
  1289. -     Other *_free_list() functions are similar.
  1290. -*/
  1291. -
  1292. -int zv_free_vars(ZVEC **pv,...)
  1293. -{
  1294. -   va_list ap;
  1295. -   int i=1;
  1296. -   ZVEC **par;
  1297. -   
  1298. -   zv_free(*pv);
  1299. -   *pv = ZVNULL;
  1300. -   va_start(ap, pv);
  1301. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1302. -      zv_free(*par); 
  1303. -      *par = ZVNULL;
  1304. -      i++;
  1305. -   } 
  1306. -
  1307. -   va_end(ap);
  1308. -   return i;
  1309. -}
  1310. -
  1311. -
  1312. -
  1313. -int zm_free_vars(ZMAT **va,...)
  1314. -{
  1315. -   va_list ap;
  1316. -   int i=1;
  1317. -   ZMAT **par;
  1318. -   
  1319. -   zm_free(*va);
  1320. -   *va = ZMNULL;
  1321. -   va_start(ap, va);
  1322. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1323. -      zm_free(*par); 
  1324. -      *par = ZMNULL;
  1325. -      i++;
  1326. -   } 
  1327. -
  1328. -   va_end(ap);
  1329. -   return i;
  1330. -}
  1331. -
  1332. -
  1333. -
  1334. -#elif VARARGS
  1335. -
  1336. -#include <varargs.h>
  1337. -
  1338. -/* To allocate memory to many arguments. 
  1339. -   The function should be called:
  1340. -   v_get_vars(dim,&x,&y,&z,...,NULL);
  1341. -   where 
  1342. -     int dim;
  1343. -     ZVEC *x, *y, *z,...;
  1344. -     The last argument should be NULL ! 
  1345. -     dim is the length of vectors x,y,z,...
  1346. -     returned value is equal to the number of allocated variables
  1347. -     Other gec_... functions are similar.
  1348. -*/
  1349. -
  1350. -int zv_get_vars(va_alist) va_dcl
  1351. -{
  1352. -   va_list ap;
  1353. -   int dim,i=0;
  1354. -   ZVEC **par;
  1355. -   
  1356. -   va_start(ap);
  1357. -   dim = va_arg(ap,int);
  1358. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1359. -      *par = zv_get(dim);
  1360. -      i++;
  1361. -   } 
  1362. -
  1363. -   va_end(ap);
  1364. -   return i;
  1365. -}
  1366. -
  1367. -
  1368. -
  1369. -int zm_get_vars(va_alist) va_dcl
  1370. -{
  1371. -   va_list ap;
  1372. -   int i=0, n, m;
  1373. -   ZMAT **par;
  1374. -   
  1375. -   va_start(ap);
  1376. -   m = va_arg(ap,int);
  1377. -   n = va_arg(ap,int);
  1378. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1379. -      *par = zm_get(m,n);
  1380. -      i++;
  1381. -   } 
  1382. -
  1383. -   va_end(ap);
  1384. -   return i;
  1385. -}
  1386. -
  1387. -
  1388. -
  1389. -/* To resize memory for many arguments. 
  1390. -   The function should be called:
  1391. -   v_resize_vars(new_dim,&x,&y,&z,...,NULL);
  1392. -   where 
  1393. -     int new_dim;
  1394. -     ZVEC *x, *y, *z,...;
  1395. -     The last argument should be NULL ! 
  1396. -     rdim is the resized length of vectors x,y,z,...
  1397. -     returned value is equal to the number of allocated variables.
  1398. -     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  1399. -     argument. 
  1400. -     Other *_resize_list() functions are similar.
  1401. -*/
  1402. -
  1403. -int zv_resize_vars(va_alist) va_dcl
  1404. -{
  1405. -   va_list ap;
  1406. -   int i=0, new_dim;
  1407. -   ZVEC **par;
  1408. -   
  1409. -   va_start(ap);
  1410. -   new_dim = va_arg(ap,int);
  1411. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1412. -      *par = zv_resize(*par,new_dim);
  1413. -      i++;
  1414. -   } 
  1415. -
  1416. -   va_end(ap);
  1417. -   return i;
  1418. -}
  1419. -
  1420. -
  1421. -int zm_resize_vars(va_alist) va_dcl
  1422. -{
  1423. -   va_list ap;
  1424. -   int i=0, m, n;
  1425. -   ZMAT **par;
  1426. -   
  1427. -   va_start(ap);
  1428. -   m = va_arg(ap,int);
  1429. -   n = va_arg(ap,int);
  1430. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1431. -      *par = zm_resize(*par,m,n);
  1432. -      i++;
  1433. -   } 
  1434. -
  1435. -   va_end(ap);
  1436. -   return i;
  1437. -}
  1438. -
  1439. -
  1440. -
  1441. -/* To deallocate memory for many arguments. 
  1442. -   The function should be called:
  1443. -   v_free_vars(&x,&y,&z,...,NULL);
  1444. -   where 
  1445. -     ZVEC *x, *y, *z,...;
  1446. -     The last argument should be NULL ! 
  1447. -     There must be at least one not NULL argument.
  1448. -     returned value is equal to the number of allocated variables.
  1449. -     Returned value of x,y,z,.. is VNULL.
  1450. -     Other *_free_list() functions are similar.
  1451. -*/
  1452. -
  1453. -int zv_free_vars(va_alist) va_dcl
  1454. -{
  1455. -   va_list ap;
  1456. -   int i=0;
  1457. -   ZVEC **par;
  1458. -   
  1459. -   va_start(ap);
  1460. -   while (par = va_arg(ap,ZVEC **)) {   /* NULL ends the list*/
  1461. -      zv_free(*par); 
  1462. -      *par = ZVNULL;
  1463. -      i++;
  1464. -   } 
  1465. -
  1466. -   va_end(ap);
  1467. -   return i;
  1468. -}
  1469. -
  1470. -
  1471. -
  1472. -int zm_free_vars(va_alist) va_dcl
  1473. -{
  1474. -   va_list ap;
  1475. -   int i=0;
  1476. -   ZMAT **par;
  1477. -   
  1478. -   va_start(ap);
  1479. -   while (par = va_arg(ap,ZMAT **)) {   /* NULL ends the list*/
  1480. -      zm_free(*par); 
  1481. -      *par = ZMNULL;
  1482. -      i++;
  1483. -   } 
  1484. -
  1485. -   va_end(ap);
  1486. -   return i;
  1487. -}
  1488. -
  1489. -
  1490. -#endif
  1491. -
  1492. //GO.SYSIN DD zmemory.c
  1493. echo zvecop.c 1>&2
  1494. sed >zvecop.c <<'//GO.SYSIN DD zvecop.c' 's/^-//'
  1495. -
  1496. -/**************************************************************************
  1497. -**
  1498. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  1499. -**
  1500. -**                 Meschach Library
  1501. -** 
  1502. -** This Meschach Library is provided "as is" without any express 
  1503. -** or implied warranty of any kind with respect to this software. 
  1504. -** In particular the authors shall not be liable for any direct, 
  1505. -** indirect, special, incidental or consequential damages arising 
  1506. -** in any way from use of the software.
  1507. -** 
  1508. -** Everyone is granted permission to copy, modify and redistribute this
  1509. -** Meschach Library, provided:
  1510. -**  1.  All copies contain this copyright notice.
  1511. -**  2.  All modified copies shall carry a notice stating who
  1512. -**      made the last modification and the date of such modification.
  1513. -**  3.  No charge is made for this software or works derived from it.  
  1514. -**      This clause shall not be construed as constraining other software
  1515. -**      distributed on the same medium as this software, nor is a
  1516. -**      distribution fee considered a charge.
  1517. -**
  1518. -***************************************************************************/
  1519. -
  1520. -
  1521. -#include    <stdio.h>
  1522. -#include    "matrix.h"
  1523. -#include    "zmatrix.h"
  1524. -static    char    rcsid[] = "$Id: zvecop.c,v 1.2 1994/03/08 05:51:07 des Exp $";
  1525. -
  1526. -
  1527. -
  1528. -/* _zin_prod -- inner product of two vectors from i0 downwards
  1529. -    -- flag != 0 means compute sum_i a[i]*.b[i];
  1530. -    -- flag == 0 means compute sum_i a[i].b[i] */
  1531. -complex    _zin_prod(a,b,i0,flag)
  1532. -ZVEC    *a,*b;
  1533. -u_int    i0, flag;
  1534. -{
  1535. -    u_int    limit;
  1536. -
  1537. -    if ( a==ZVNULL || b==ZVNULL )
  1538. -        error(E_NULL,"_zin_prod");
  1539. -    limit = min(a->dim,b->dim);
  1540. -    if ( i0 > limit )
  1541. -        error(E_BOUNDS,"_zin_prod");
  1542. -
  1543. -    return __zip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0),flag);
  1544. -}
  1545. -
  1546. -/* zv_mlt -- scalar-vector multiply -- may be in-situ */
  1547. -ZVEC    *zv_mlt(scalar,vector,out)
  1548. -complex    scalar;
  1549. -ZVEC    *vector,*out;
  1550. -{
  1551. -    /* u_int    dim, i; */
  1552. -    /* complex    *out_ve, *vec_ve; */
  1553. -
  1554. -    if ( vector==ZVNULL )
  1555. -        error(E_NULL,"zv_mlt");
  1556. -    if ( out==ZVNULL || out->dim != vector->dim )
  1557. -        out = zv_resize(out,vector->dim);
  1558. -    if ( scalar.re == 0.0 && scalar.im == 0.0 )
  1559. -        return zv_zero(out);
  1560. -    if ( scalar.re == 1.0 && scalar.im == 0.0 )
  1561. -        return zv_copy(vector,out);
  1562. -
  1563. -    __zmlt__(vector->ve,scalar,out->ve,(int)(vector->dim));
  1564. -
  1565. -    return (out);
  1566. -}
  1567. -
  1568. -/* zv_add -- vector addition -- may be in-situ */
  1569. -ZVEC    *zv_add(vec1,vec2,out)
  1570. -ZVEC    *vec1,*vec2,*out;
  1571. -{
  1572. -    u_int    dim;
  1573. -
  1574. -    if ( vec1==ZVNULL || vec2==ZVNULL )
  1575. -        error(E_NULL,"zv_add");
  1576. -    if ( vec1->dim != vec2->dim )
  1577. -        error(E_SIZES,"zv_add");
  1578. -    if ( out==ZVNULL || out->dim != vec1->dim )
  1579. -        out = zv_resize(out,vec1->dim);
  1580. -    dim = vec1->dim;
  1581. -    __zadd__(vec1->ve,vec2->ve,out->ve,(int)dim);
  1582. -
  1583. -    return (out);
  1584. -}
  1585. -
  1586. -/* zv_mltadd -- scalar/vector multiplication and addition
  1587. -        -- out = v1 + scale.v2        */
  1588. -ZVEC    *zv_mltadd(v1,v2,scale,out)
  1589. -ZVEC    *v1,*v2,*out;
  1590. -complex    scale;
  1591. -{
  1592. -    /* register u_int    dim, i; */
  1593. -    /* complex    *out_ve, *v1_ve, *v2_ve; */
  1594. -
  1595. -    if ( v1==ZVNULL || v2==ZVNULL )
  1596. -        error(E_NULL,"zv_mltadd");
  1597. -    if ( v1->dim != v2->dim )
  1598. -        error(E_SIZES,"zv_mltadd");
  1599. -    if ( scale.re == 0.0 && scale.im == 0.0 )
  1600. -        return zv_copy(v1,out);
  1601. -    if ( scale.re == 1.0 && scale.im == 0.0 )
  1602. -        return zv_add(v1,v2,out);
  1603. -
  1604. -    if ( v2 != out )
  1605. -    {
  1606. -        tracecatch(out = zv_copy(v1,out),"zv_mltadd");
  1607. -
  1608. -        /* dim = v1->dim; */
  1609. -        __zmltadd__(out->ve,v2->ve,scale,(int)(v1->dim),0);
  1610. -    }
  1611. -    else
  1612. -    {
  1613. -        tracecatch(out = zv_mlt(scale,v2,out),"zv_mltadd");
  1614. -        out = zv_add(v1,out,out);
  1615. -    }
  1616. -
  1617. -    return (out);
  1618. -}
  1619. -
  1620. -/* zv_sub -- vector subtraction -- may be in-situ */
  1621. -ZVEC    *zv_sub(vec1,vec2,out)
  1622. -ZVEC    *vec1,*vec2,*out;
  1623. -{
  1624. -    /* u_int    i, dim; */
  1625. -    /* complex    *out_ve, *vec1_ve, *vec2_ve; */
  1626. -
  1627. -    if ( vec1==ZVNULL || vec2==ZVNULL )
  1628. -        error(E_NULL,"zv_sub");
  1629. -    if ( vec1->dim != vec2->dim )
  1630. -        error(E_SIZES,"zv_sub");
  1631. -    if ( out==ZVNULL || out->dim != vec1->dim )
  1632. -        out = zv_resize(out,vec1->dim);
  1633. -
  1634. -    __zsub__(vec1->ve,vec2->ve,out->ve,(int)(vec1->dim));
  1635. -
  1636. -    return (out);
  1637. -}
  1638. -
  1639. -/* zv_map -- maps function f over components of x: out[i] = f(x[i])
  1640. -    -- _zv_map sets out[i] = f(x[i],params) */
  1641. -ZVEC    *zv_map(f,x,out)
  1642. -#ifdef PROTOYPES_IN_STRUCT
  1643. -complex    (*f)(complex);
  1644. -#else
  1645. -complex (*f)();
  1646. -#endif
  1647. -ZVEC    *x, *out;
  1648. -{
  1649. -    complex    *x_ve, *out_ve;
  1650. -    int    i, dim;
  1651. -
  1652. -    if ( ! x || ! f )
  1653. -        error(E_NULL,"zv_map");
  1654. -    if ( ! out || out->dim != x->dim )
  1655. -        out = zv_resize(out,x->dim);
  1656. -
  1657. -    dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  1658. -    for ( i = 0; i < dim; i++ )
  1659. -        out_ve[i] = (*f)(x_ve[i]);
  1660. -
  1661. -    return out;
  1662. -}
  1663. -
  1664. -ZVEC    *_zv_map(f,params,x,out)
  1665. -#ifdef PROTOTYPES_IN_STRUCT
  1666. -complex    (*f)(void *,complex);
  1667. -#else
  1668. -complex    (*f)();
  1669. -#endif
  1670. -ZVEC    *x, *out;
  1671. -void    *params;
  1672. -{
  1673. -    complex    *x_ve, *out_ve;
  1674. -    int    i, dim;
  1675. -
  1676. -    if ( ! x || ! f )
  1677. -        error(E_NULL,"_zv_map");
  1678. -    if ( ! out || out->dim != x->dim )
  1679. -        out = zv_resize(out,x->dim);
  1680. -
  1681. -    dim = x->dim;    x_ve = x->ve;    out_ve = out->ve;
  1682. -    for ( i = 0; i < dim; i++ )
  1683. -        out_ve[i] = (*f)(params,x_ve[i]);
  1684. -
  1685. -    return out;
  1686. -}
  1687. -
  1688. -/* zv_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */
  1689. -ZVEC    *zv_lincomb(n,v,a,out)
  1690. -int    n;    /* number of a's and v's */
  1691. -complex    a[];
  1692. -ZVEC    *v[], *out;
  1693. -{
  1694. -    int    i;
  1695. -
  1696. -    if ( ! a || ! v )
  1697. -        error(E_NULL,"zv_lincomb");
  1698. -    if ( n <= 0 )
  1699. -        return ZVNULL;
  1700. -
  1701. -    for ( i = 1; i < n; i++ )
  1702. -        if ( out == v[i] )
  1703. -            error(E_INSITU,"zv_lincomb");
  1704. -
  1705. -    out = zv_mlt(a[0],v[0],out);
  1706. -    for ( i = 1; i < n; i++ )
  1707. -    {
  1708. -        if ( ! v[i] )
  1709. -            error(E_NULL,"zv_lincomb");
  1710. -        if ( v[i]->dim != out->dim )
  1711. -            error(E_SIZES,"zv_lincomb");
  1712. -        out = zv_mltadd(out,v[i],a[i],out);
  1713. -    }
  1714. -
  1715. -    return out;
  1716. -}
  1717. -
  1718. -
  1719. -#ifdef ANSI_C
  1720. -
  1721. -
  1722. -/* zv_linlist -- linear combinations taken from a list of arguments;
  1723. -   calling:
  1724. -      zv_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  1725. -   where vi are vectors (ZVEC *) and ai are numbers (complex)
  1726. -*/
  1727. -
  1728. -ZVEC    *zv_linlist(ZVEC *out,ZVEC *v1,complex a1,...)
  1729. -{
  1730. -   va_list ap;
  1731. -   ZVEC *par;
  1732. -   complex a_par;
  1733. -
  1734. -   if ( ! v1 )
  1735. -     return ZVNULL;
  1736. -   
  1737. -   va_start(ap, a1);
  1738. -   out = zv_mlt(a1,v1,out);
  1739. -   
  1740. -   while (par = va_arg(ap,ZVEC *)) {   /* NULL ends the list*/
  1741. -      a_par = va_arg(ap,complex);
  1742. -      if (a_par.re == 0.0 && a_par.im == 0.0) continue;
  1743. -      if ( out == par )        
  1744. -    error(E_INSITU,"zv_linlist");
  1745. -      if ( out->dim != par->dim )    
  1746. -    error(E_SIZES,"zv_linlist");
  1747. -
  1748. -      if (a_par.re == 1.0 && a_par.im == 0.0)
  1749. -    out = zv_add(out,par,out);
  1750. -      else if (a_par.re == -1.0 && a_par.im == 0.0)
  1751. -    out = zv_sub(out,par,out);
  1752. -      else
  1753. -    out = zv_mltadd(out,par,a_par,out); 
  1754. -   } 
  1755. -   
  1756. -   va_end(ap);
  1757. -   return out;
  1758. -}
  1759. -
  1760. -
  1761. -#elif VARARGS
  1762. -
  1763. -/* zv_linlist -- linear combinations taken from a list of arguments;
  1764. -   calling:
  1765. -      zv_linlist(out,v1,a1,v2,a2,...,vn,an,NULL);
  1766. -   where vi are vectors (ZVEC *) and ai are numbers (complex)
  1767. -*/
  1768. -ZVEC  *zv_linlist(va_alist) va_dcl
  1769. -{
  1770. -   va_list ap;
  1771. -   ZVEC *par, *out;
  1772. -   complex a_par;
  1773. -
  1774. -   va_start(ap);
  1775. -   out = va_arg(ap,ZVEC *);
  1776. -   par = va_arg(ap,ZVEC *);
  1777. -   if ( ! par ) {
  1778. -      va_end(ap);
  1779. -      return ZVNULL;
  1780. -   }
  1781. -   
  1782. -   a_par = va_arg(ap,complex);
  1783. -   out = zv_mlt(a_par,par,out);
  1784. -   
  1785. -   while (par = va_arg(ap,ZVEC *)) {   /* NULL ends the list*/
  1786. -      a_par = va_arg(ap,complex);
  1787. -      if (a_par.re == 0.0 && a_par.im == 0.0) continue;
  1788. -      if ( out == par )        
  1789. -    error(E_INSITU,"zv_linlist");
  1790. -      if ( out->dim != par->dim )    
  1791. -    error(E_SIZES,"zv_linlist");
  1792. -
  1793. -      if (a_par.re == 1.0 && a_par.im == 0.0)
  1794. -    out = zv_add(out,par,out);
  1795. -      else if (a_par.re == -1.0 && a_par.im == 0.0)
  1796. -    out = zv_sub(out,par,out);
  1797. -      else
  1798. -    out = zv_mltadd(out,par,a_par,out); 
  1799. -   } 
  1800. -   
  1801. -   va_end(ap);
  1802. -   return out;
  1803. -}
  1804. -
  1805. -
  1806. -#endif
  1807. -
  1808. -
  1809. -
  1810. -/* zv_star -- computes componentwise (Hadamard) product of x1 and x2
  1811. -    -- result out is returned */
  1812. -ZVEC    *zv_star(x1, x2, out)
  1813. -ZVEC    *x1, *x2, *out;
  1814. -{
  1815. -    int        i;
  1816. -    Real    t_re, t_im;
  1817. -
  1818. -    if ( ! x1 || ! x2 )
  1819. -    error(E_NULL,"zv_star");
  1820. -    if ( x1->dim != x2->dim )
  1821. -    error(E_SIZES,"zv_star");
  1822. -    out = zv_resize(out,x1->dim);
  1823. -
  1824. -    for ( i = 0; i < x1->dim; i++ )
  1825. -    {
  1826. -    /* out->ve[i] = x1->ve[i] * x2->ve[i]; */
  1827. -    t_re = x1->ve[i].re*x2->ve[i].re - x1->ve[i].im*x2->ve[i].im;
  1828. -    t_im = x1->ve[i].re*x2->ve[i].im + x1->ve[i].im*x2->ve[i].re;
  1829. -    out->ve[i].re = t_re;
  1830. -    out->ve[i].im = t_im;
  1831. -    }
  1832. -
  1833. -    return out;
  1834. -}
  1835. -
  1836. -/* zv_slash -- computes componentwise ratio of x2 and x1
  1837. -    -- out[i] = x2[i] / x1[i]
  1838. -    -- if x1[i] == 0 for some i, then raise E_SING error
  1839. -    -- result out is returned */
  1840. -ZVEC    *zv_slash(x1, x2, out)
  1841. -ZVEC    *x1, *x2, *out;
  1842. -{
  1843. -    int        i;
  1844. -    Real    r2, t_re, t_im;
  1845. -    complex    tmp;
  1846. -
  1847. -    if ( ! x1 || ! x2 )
  1848. -    error(E_NULL,"zv_slash");
  1849. -    if ( x1->dim != x2->dim )
  1850. -    error(E_SIZES,"zv_slash");
  1851. -    out = zv_resize(out,x1->dim);
  1852. -
  1853. -    for ( i = 0; i < x1->dim; i++ )
  1854. -    {
  1855. -    r2 = x1->ve[i].re*x1->ve[i].re + x1->ve[i].im*x1->ve[i].im;
  1856. -    if ( r2 == 0.0 )
  1857. -        error(E_SING,"zv_slash");
  1858. -    tmp.re =   x1->ve[i].re / r2;
  1859. -    tmp.im = - x1->ve[i].im / r2;
  1860. -    t_re = tmp.re*x2->ve[i].re - tmp.im*x2->ve[i].im;
  1861. -    t_im = tmp.re*x2->ve[i].im - tmp.im*x2->ve[i].re;
  1862. -    out->ve[i].re = t_re;
  1863. -    out->ve[i].im = t_im;
  1864. -    }
  1865. -
  1866. -    return out;
  1867. -}
  1868. -
  1869. -/* zv_sum -- returns sum of entries of a vector */
  1870. -complex    zv_sum(x)
  1871. -ZVEC    *x;
  1872. -{
  1873. -    int        i;
  1874. -    complex    sum;
  1875. -
  1876. -    if ( ! x )
  1877. -    error(E_NULL,"zv_sum");
  1878. -
  1879. -    sum.re = sum.im = 0.0;
  1880. -    for ( i = 0; i < x->dim; i++ )
  1881. -    {
  1882. -    sum.re += x->ve[i].re;
  1883. -    sum.im += x->ve[i].im;
  1884. -    }
  1885. -
  1886. -    return sum;
  1887. -}
  1888. -
  1889. -/* px_zvec -- permute vector */
  1890. -ZVEC    *px_zvec(px,vector,out)
  1891. -PERM    *px;
  1892. -ZVEC    *vector,*out;
  1893. -{
  1894. -    u_int    old_i, i, size, start;
  1895. -    complex    tmp;
  1896. -    
  1897. -    if ( px==PNULL || vector==ZVNULL )
  1898. -    error(E_NULL,"px_zvec");
  1899. -    if ( px->size > vector->dim )
  1900. -    error(E_SIZES,"px_zvec");
  1901. -    if ( out==ZVNULL || out->dim < vector->dim )
  1902. -    out = zv_resize(out,vector->dim);
  1903. -    
  1904. -    size = px->size;
  1905. -    if ( size == 0 )
  1906. -    return zv_copy(vector,out);
  1907. -    
  1908. -    if ( out != vector )
  1909. -    {
  1910. -    for ( i=0; i<size; i++ )
  1911. -        if ( px->pe[i] >= size )
  1912. -        error(E_BOUNDS,"px_vec");
  1913. -        else
  1914. -        out->ve[i] = vector->ve[px->pe[i]];
  1915. -    }
  1916. -    else
  1917. -    {    /* in situ algorithm */
  1918. -    start = 0;
  1919. -    while ( start < size )
  1920. -    {
  1921. -        old_i = start;
  1922. -        i = px->pe[old_i];
  1923. -        if ( i >= size )
  1924. -        {
  1925. -        start++;
  1926. -        continue;
  1927. -        }
  1928. -        tmp = vector->ve[start];
  1929. -        while ( TRUE )
  1930. -        {
  1931. -        vector->ve[old_i] = vector->ve[i];
  1932. -        px->pe[old_i] = i+size;
  1933. -        old_i = i;
  1934. -        i = px->pe[old_i];
  1935. -        if ( i >= size )
  1936. -            break;
  1937. -        if ( i == start )
  1938. -        {
  1939. -            vector->ve[old_i] = tmp;
  1940. -            px->pe[old_i] = i+size;
  1941. -            break;
  1942. -        }
  1943. -        }
  1944. -        start++;
  1945. -    }
  1946. -    
  1947. -    for ( i = 0; i < size; i++ )
  1948. -        if ( px->pe[i] < size )
  1949. -        error(E_BOUNDS,"px_vec");
  1950. -        else
  1951. -        px->pe[i] = px->pe[i]-size;
  1952. -    }
  1953. -    
  1954. -    return out;
  1955. -}
  1956. -
  1957. -/* pxinv_zvec -- apply the inverse of px to x, returning the result in out
  1958. -        -- may NOT be in situ */
  1959. -ZVEC    *pxinv_zvec(px,x,out)
  1960. -PERM    *px;
  1961. -ZVEC    *x, *out;
  1962. -{
  1963. -    u_int    i, size;
  1964. -    
  1965. -    if ( ! px || ! x )
  1966. -    error(E_NULL,"pxinv_zvec");
  1967. -    if ( px->size > x->dim )
  1968. -    error(E_SIZES,"pxinv_zvec");
  1969. -    if ( ! out || out->dim < x->dim )
  1970. -    out = zv_resize(out,x->dim);
  1971. -    
  1972. -    size = px->size;
  1973. -    if ( size == 0 )
  1974. -    return zv_copy(x,out);
  1975. -    if ( out != x )
  1976. -    {
  1977. -    for ( i=0; i<size; i++ )
  1978. -        if ( px->pe[i] >= size )
  1979. -        error(E_BOUNDS,"pxinv_vec");
  1980. -        else
  1981. -        out->ve[px->pe[i]] = x->ve[i];
  1982. -    }
  1983. -    else
  1984. -    {    /* in situ algorithm --- cheat's way out */
  1985. -    px_inv(px,px);
  1986. -    px_zvec(px,x,out);
  1987. -    px_inv(px,px);
  1988. -    }
  1989. -    
  1990. -    
  1991. -    return out;
  1992. -}
  1993. -
  1994. -/* zv_rand -- randomise a complex vector; uniform in [0,1)+[0,1)*i */
  1995. -ZVEC    *zv_rand(x)
  1996. -ZVEC    *x;
  1997. -{
  1998. -    if ( ! x )
  1999. -    error(E_NULL,"zv_rand");
  2000. -
  2001. -    mrandlist((Real *)(x->ve),2*x->dim);
  2002. -
  2003. -    return x;
  2004. -}
  2005. //GO.SYSIN DD zvecop.c
  2006. echo zmatop.c 1>&2
  2007. sed >zmatop.c <<'//GO.SYSIN DD zmatop.c' 's/^-//'
  2008. -
  2009. -/**************************************************************************
  2010. -**
  2011. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2012. -**
  2013. -**                 Meschach Library
  2014. -** 
  2015. -** This Meschach Library is provided "as is" without any express 
  2016. -** or implied warranty of any kind with respect to this software. 
  2017. -** In particular the authors shall not be liable for any direct, 
  2018. -** indirect, special, incidental or consequential damages arising 
  2019. -** in any way from use of the software.
  2020. -** 
  2021. -** Everyone is granted permission to copy, modify and redistribute this
  2022. -** Meschach Library, provided:
  2023. -**  1.  All copies contain this copyright notice.
  2024. -**  2.  All modified copies shall carry a notice stating who
  2025. -**      made the last modification and the date of such modification.
  2026. -**  3.  No charge is made for this software or works derived from it.  
  2027. -**      This clause shall not be construed as constraining other software
  2028. -**      distributed on the same medium as this software, nor is a
  2029. -**      distribution fee considered a charge.
  2030. -**
  2031. -***************************************************************************/
  2032. -
  2033. -
  2034. -
  2035. -#include    <stdio.h>
  2036. -#include    "zmatrix.h"
  2037. -
  2038. -static    char    rcsid[] = "$Id: zmatop.c,v 1.1 1994/01/13 04:24:46 des Exp $";
  2039. -
  2040. -
  2041. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  2042. -
  2043. -/* zm_add -- matrix addition -- may be in-situ */
  2044. -ZMAT    *zm_add(mat1,mat2,out)
  2045. -ZMAT    *mat1,*mat2,*out;
  2046. -{
  2047. -    u_int    m,n,i;
  2048. -    
  2049. -    if ( mat1==ZMNULL || mat2==ZMNULL )
  2050. -    error(E_NULL,"zm_add");
  2051. -    if ( mat1->m != mat2->m || mat1->n != mat2->n )
  2052. -    error(E_SIZES,"zm_add");
  2053. -    if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  2054. -    out = zm_resize(out,mat1->m,mat1->n);
  2055. -    m = mat1->m;    n = mat1->n;
  2056. -    for ( i=0; i<m; i++ )
  2057. -    {
  2058. -    __zadd__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  2059. -    /**************************************************
  2060. -      for ( j=0; j<n; j++ )
  2061. -      out->me[i][j] = mat1->me[i][j]+mat2->me[i][j];
  2062. -      **************************************************/
  2063. -    }
  2064. -    
  2065. -    return (out);
  2066. -}
  2067. -
  2068. -/* zm_sub -- matrix subtraction -- may be in-situ */
  2069. -ZMAT    *zm_sub(mat1,mat2,out)
  2070. -ZMAT    *mat1,*mat2,*out;
  2071. -{
  2072. -    u_int    m,n,i;
  2073. -    
  2074. -    if ( mat1==ZMNULL || mat2==ZMNULL )
  2075. -    error(E_NULL,"zm_sub");
  2076. -    if ( mat1->m != mat2->m || mat1->n != mat2->n )
  2077. -    error(E_SIZES,"zm_sub");
  2078. -    if ( out==ZMNULL || out->m != mat1->m || out->n != mat1->n )
  2079. -    out = zm_resize(out,mat1->m,mat1->n);
  2080. -    m = mat1->m;    n = mat1->n;
  2081. -    for ( i=0; i<m; i++ )
  2082. -    {
  2083. -    __zsub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
  2084. -    /**************************************************
  2085. -      for ( j=0; j<n; j++ )
  2086. -      out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
  2087. -    **************************************************/
  2088. -    }
  2089. -    
  2090. -    return (out);
  2091. -}
  2092. -
  2093. -/*
  2094. -  Note: In the following routines, "adjoint" means complex conjugate
  2095. -  transpose:
  2096. -  A* = conjugate(A^T)
  2097. -  */
  2098. -
  2099. -/* zm_mlt -- matrix-matrix multiplication */
  2100. -ZMAT    *zm_mlt(A,B,OUT)
  2101. -ZMAT    *A,*B,*OUT;
  2102. -{
  2103. -    u_int    i, /* j, */ k, m, n, p;
  2104. -    complex    **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
  2105. -    
  2106. -    if ( A==ZMNULL || B==ZMNULL )
  2107. -    error(E_NULL,"zm_mlt");
  2108. -    if ( A->n != B->m )
  2109. -    error(E_SIZES,"zm_mlt");
  2110. -    if ( A == OUT || B == OUT )
  2111. -    error(E_INSITU,"zm_mlt");
  2112. -    m = A->m;    n = A->n;    p = B->n;
  2113. -    A_v = A->me;        B_v = B->me;
  2114. -    
  2115. -    if ( OUT==ZMNULL || OUT->m != A->m || OUT->n != B->n )
  2116. -    OUT = zm_resize(OUT,A->m,B->n);
  2117. -    
  2118. -    /****************************************************************
  2119. -      for ( i=0; i<m; i++ )
  2120. -      for  ( j=0; j<p; j++ )
  2121. -      {
  2122. -      sum = 0.0;
  2123. -      for ( k=0; k<n; k++ )
  2124. -      sum += A_v[i][k]*B_v[k][j];
  2125. -      OUT->me[i][j] = sum;
  2126. -      }
  2127. -    ****************************************************************/
  2128. -    zm_zero(OUT);
  2129. -    for ( i=0; i<m; i++ )
  2130. -    for ( k=0; k<n; k++ )
  2131. -    {
  2132. -        if ( ! is_zero(A_v[i][k]) )
  2133. -        __zmltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p,Z_NOCONJ);
  2134. -        /**************************************************
  2135. -          B_row = B_v[k];    OUT_row = OUT->me[i];
  2136. -          for ( j=0; j<p; j++ )
  2137. -          (*OUT_row++) += tmp*(*B_row++);
  2138. -        **************************************************/
  2139. -    }
  2140. -    
  2141. -    return OUT;
  2142. -}
  2143. -
  2144. -/* zmma_mlt -- matrix-matrix adjoint multiplication
  2145. -   -- A.B* is returned, and stored in OUT */
  2146. -ZMAT    *zmma_mlt(A,B,OUT)
  2147. -ZMAT    *A, *B, *OUT;
  2148. -{
  2149. -    int    i, j, limit;
  2150. -    /* complex    *A_row, *B_row, sum; */
  2151. -    
  2152. -    if ( ! A || ! B )
  2153. -    error(E_NULL,"zmma_mlt");
  2154. -    if ( A == OUT || B == OUT )
  2155. -    error(E_INSITU,"zmma_mlt");
  2156. -    if ( A->n != B->n )
  2157. -    error(E_SIZES,"zmma_mlt");
  2158. -    if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
  2159. -    OUT = zm_resize(OUT,A->m,B->m);
  2160. -    
  2161. -    limit = A->n;
  2162. -    for ( i = 0; i < A->m; i++ )
  2163. -    for ( j = 0; j < B->m; j++ )
  2164. -    {
  2165. -        OUT->me[i][j] = __zip__(B->me[j],A->me[i],(int)limit,Z_CONJ);
  2166. -        /**************************************************
  2167. -          sum = 0.0;
  2168. -          A_row = A->me[i];
  2169. -          B_row = B->me[j];
  2170. -          for ( k = 0; k < limit; k++ )
  2171. -          sum += (*A_row++)*(*B_row++);
  2172. -          OUT->me[i][j] = sum;
  2173. -          **************************************************/
  2174. -    }
  2175. -    
  2176. -    return OUT;
  2177. -}
  2178. -
  2179. -/* zmam_mlt -- matrix adjoint-matrix multiplication
  2180. -   -- A*.B is returned, result stored in OUT */
  2181. -ZMAT    *zmam_mlt(A,B,OUT)
  2182. -ZMAT    *A, *B, *OUT;
  2183. -{
  2184. -    int    i, k, limit;
  2185. -    /* complex    *B_row, *OUT_row, multiplier; */
  2186. -    complex    tmp;
  2187. -    
  2188. -    if ( ! A || ! B )
  2189. -    error(E_NULL,"zmam_mlt");
  2190. -    if ( A == OUT || B == OUT )
  2191. -    error(E_INSITU,"zmam_mlt");
  2192. -    if ( A->m != B->m )
  2193. -    error(E_SIZES,"zmam_mlt");
  2194. -    if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
  2195. -    OUT = zm_resize(OUT,A->n,B->n);
  2196. -    
  2197. -    limit = B->n;
  2198. -    zm_zero(OUT);
  2199. -    for ( k = 0; k < A->m; k++ )
  2200. -    for ( i = 0; i < A->n; i++ )
  2201. -    {
  2202. -        tmp.re =   A->me[k][i].re;
  2203. -        tmp.im = - A->me[k][i].im;
  2204. -        if ( ! is_zero(tmp) )
  2205. -        __zmltadd__(OUT->me[i],B->me[k],tmp,(int)limit,Z_NOCONJ);
  2206. -    }
  2207. -    
  2208. -    return OUT;
  2209. -}
  2210. -
  2211. -/* zmv_mlt -- matrix-vector multiplication 
  2212. -   -- Note: b is treated as a column vector */
  2213. -ZVEC    *zmv_mlt(A,b,out)
  2214. -ZMAT    *A;
  2215. -ZVEC    *b,*out;
  2216. -{
  2217. -    u_int    i, m, n;
  2218. -    complex    **A_v, *b_v /*, *A_row */;
  2219. -    /* register complex    sum; */
  2220. -    
  2221. -    if ( A==ZMNULL || b==ZVNULL )
  2222. -    error(E_NULL,"zmv_mlt");
  2223. -    if ( A->n != b->dim )
  2224. -    error(E_SIZES,"zmv_mlt");
  2225. -    if ( b == out )
  2226. -    error(E_INSITU,"zmv_mlt");
  2227. -    if ( out == ZVNULL || out->dim != A->m )
  2228. -    out = zv_resize(out,A->m);
  2229. -    
  2230. -    m = A->m;        n = A->n;
  2231. -    A_v = A->me;        b_v = b->ve;
  2232. -    for ( i=0; i<m; i++ )
  2233. -    {
  2234. -    /* for ( j=0; j<n; j++ )
  2235. -       sum += A_v[i][j]*b_v[j]; */
  2236. -    out->ve[i] = __zip__(A_v[i],b_v,(int)n,Z_NOCONJ);
  2237. -    /**************************************************
  2238. -      A_row = A_v[i];        b_v = b->ve;
  2239. -      for ( j=0; j<n; j++ )
  2240. -      sum += (*A_row++)*(*b_v++);
  2241. -      out->ve[i] = sum;
  2242. -    **************************************************/
  2243. -    }
  2244. -    
  2245. -    return out;
  2246. -}
  2247. -
  2248. -/* zsm_mlt -- scalar-matrix multiply -- may be in-situ */
  2249. -ZMAT    *zsm_mlt(scalar,matrix,out)
  2250. -complex    scalar;
  2251. -ZMAT    *matrix,*out;
  2252. -{
  2253. -    u_int    m,n,i;
  2254. -    
  2255. -    if ( matrix==ZMNULL )
  2256. -    error(E_NULL,"zsm_mlt");
  2257. -    if ( out==ZMNULL || out->m != matrix->m || out->n != matrix->n )
  2258. -    out = zm_resize(out,matrix->m,matrix->n);
  2259. -    m = matrix->m;    n = matrix->n;
  2260. -    for ( i=0; i<m; i++ )
  2261. -    __zmlt__(matrix->me[i],scalar,out->me[i],(int)n);
  2262. -    /**************************************************
  2263. -      for ( j=0; j<n; j++ )
  2264. -      out->me[i][j] = scalar*matrix->me[i][j];
  2265. -      **************************************************/
  2266. -    return (out);
  2267. -}
  2268. -
  2269. -/* zvm_mlt -- vector adjoint-matrix multiplication */
  2270. -ZVEC    *zvm_mlt(A,b,out)
  2271. -ZMAT    *A;
  2272. -ZVEC    *b,*out;
  2273. -{
  2274. -    u_int    j,m,n;
  2275. -    /* complex    sum,**A_v,*b_v; */
  2276. -    
  2277. -    if ( A==ZMNULL || b==ZVNULL )
  2278. -    error(E_NULL,"zvm_mlt");
  2279. -    if ( A->m != b->dim )
  2280. -    error(E_SIZES,"zvm_mlt");
  2281. -    if ( b == out )
  2282. -    error(E_INSITU,"zvm_mlt");
  2283. -    if ( out == ZVNULL || out->dim != A->n )
  2284. -    out = zv_resize(out,A->n);
  2285. -    
  2286. -    m = A->m;        n = A->n;
  2287. -    
  2288. -    zv_zero(out);
  2289. -    for ( j = 0; j < m; j++ )
  2290. -    if ( b->ve[j].re != 0.0 || b->ve[j].im != 0.0  )
  2291. -        __zmltadd__(out->ve,A->me[j],b->ve[j],(int)n,Z_CONJ);
  2292. -    /**************************************************
  2293. -      A_v = A->me;        b_v = b->ve;
  2294. -      for ( j=0; j<n; j++ )
  2295. -      {
  2296. -      sum = 0.0;
  2297. -      for ( i=0; i<m; i++ )
  2298. -      sum += b_v[i]*A_v[i][j];
  2299. -      out->ve[j] = sum;
  2300. -      }
  2301. -      **************************************************/
  2302. -    
  2303. -    return out;
  2304. -}
  2305. -
  2306. -/* zm_adjoint -- adjoint matrix */
  2307. -ZMAT    *zm_adjoint(in,out)
  2308. -ZMAT    *in, *out;
  2309. -{
  2310. -    int    i, j;
  2311. -    int    in_situ;
  2312. -    complex    tmp;
  2313. -    
  2314. -    if ( in == ZMNULL )
  2315. -    error(E_NULL,"zm_adjoint");
  2316. -    if ( in == out && in->n != in->m )
  2317. -    error(E_INSITU2,"zm_adjoint");
  2318. -    in_situ = ( in == out );
  2319. -    if ( out == ZMNULL || out->m != in->n || out->n != in->m )
  2320. -    out = zm_resize(out,in->n,in->m);
  2321. -    
  2322. -    if ( ! in_situ )
  2323. -    {
  2324. -    for ( i = 0; i < in->m; i++ )
  2325. -        for ( j = 0; j < in->n; j++ )
  2326. -        {
  2327. -        out->me[j][i].re =   in->me[i][j].re;
  2328. -        out->me[j][i].im = - in->me[i][j].im;
  2329. -        }
  2330. -    }
  2331. -    else
  2332. -    {
  2333. -    for ( i = 0 ; i < in->m; i++ )
  2334. -    {
  2335. -        for ( j = 0; j < i; j++ )
  2336. -        {
  2337. -        tmp.re = in->me[i][j].re;
  2338. -        tmp.im = in->me[i][j].im;
  2339. -        in->me[i][j].re =   in->me[j][i].re;
  2340. -        in->me[i][j].im = - in->me[j][i].im;
  2341. -        in->me[j][i].re =   tmp.re;
  2342. -        in->me[j][i].im = - tmp.im;
  2343. -        }
  2344. -        in->me[i][i].im = - in->me[i][i].im;
  2345. -    }
  2346. -    }
  2347. -    
  2348. -    return out;
  2349. -}
  2350. -
  2351. -/* zswap_rows -- swaps rows i and j of matrix A upto column lim */
  2352. -ZMAT    *zswap_rows(A,i,j,lo,hi)
  2353. -ZMAT    *A;
  2354. -int    i, j, lo, hi;
  2355. -{
  2356. -    int    k;
  2357. -    complex    **A_me, tmp;
  2358. -    
  2359. -    if ( ! A )
  2360. -    error(E_NULL,"swap_rows");
  2361. -    if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
  2362. -    error(E_SIZES,"swap_rows");
  2363. -    lo = max(0,lo);
  2364. -    hi = min(hi,A->n-1);
  2365. -    A_me = A->me;
  2366. -    
  2367. -    for ( k = lo; k <= hi; k++ )
  2368. -    {
  2369. -    tmp = A_me[k][i];
  2370. -    A_me[k][i] = A_me[k][j];
  2371. -    A_me[k][j] = tmp;
  2372. -    }
  2373. -    return A;
  2374. -}
  2375. -
  2376. -/* zswap_cols -- swap columns i and j of matrix A upto row lim */
  2377. -ZMAT    *zswap_cols(A,i,j,lo,hi)
  2378. -ZMAT    *A;
  2379. -int    i, j, lo, hi;
  2380. -{
  2381. -    int    k;
  2382. -    complex    **A_me, tmp;
  2383. -    
  2384. -    if ( ! A )
  2385. -    error(E_NULL,"swap_cols");
  2386. -    if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
  2387. -    error(E_SIZES,"swap_cols");
  2388. -    lo = max(0,lo);
  2389. -    hi = min(hi,A->m-1);
  2390. -    A_me = A->me;
  2391. -    
  2392. -    for ( k = lo; k <= hi; k++ )
  2393. -    {
  2394. -    tmp = A_me[i][k];
  2395. -    A_me[i][k] = A_me[j][k];
  2396. -    A_me[j][k] = tmp;
  2397. -    }
  2398. -    return A;
  2399. -}
  2400. -
  2401. -/* mz_mltadd -- matrix-scalar multiply and add
  2402. -   -- may be in situ
  2403. -   -- returns out == A1 + s*A2 */
  2404. -ZMAT    *mz_mltadd(A1,A2,s,out)
  2405. -ZMAT    *A1, *A2, *out;
  2406. -complex    s;
  2407. -{
  2408. -    /* register complex    *A1_e, *A2_e, *out_e; */
  2409. -    /* register int    j; */
  2410. -    int    i, m, n;
  2411. -    
  2412. -    if ( ! A1 || ! A2 )
  2413. -    error(E_NULL,"mz_mltadd");
  2414. -    if ( A1->m != A2->m || A1->n != A2->n )
  2415. -    error(E_SIZES,"mz_mltadd");
  2416. -    
  2417. -    if ( s.re == 0.0 && s.im == 0.0 )
  2418. -    return zm_copy(A1,out);
  2419. -    if ( s.re == 1.0 && s.im == 0.0 )
  2420. -    return zm_add(A1,A2,out);
  2421. -    
  2422. -    tracecatch(out = zm_copy(A1,out),"mz_mltadd");
  2423. -    
  2424. -    m = A1->m;    n = A1->n;
  2425. -    for ( i = 0; i < m; i++ )
  2426. -    {
  2427. -    __zmltadd__(out->me[i],A2->me[i],s,(int)n,Z_NOCONJ);
  2428. -    /**************************************************
  2429. -      A1_e = A1->me[i];
  2430. -      A2_e = A2->me[i];
  2431. -      out_e = out->me[i];
  2432. -      for ( j = 0; j < n; j++ )
  2433. -      out_e[j] = A1_e[j] + s*A2_e[j];
  2434. -      **************************************************/
  2435. -    }
  2436. -    
  2437. -    return out;
  2438. -}
  2439. -
  2440. -/* zmv_mltadd -- matrix-vector multiply and add
  2441. -   -- may not be in situ
  2442. -   -- returns out == v1 + alpha*A*v2 */
  2443. -ZVEC    *zmv_mltadd(v1,v2,A,alpha,out)
  2444. -ZVEC    *v1, *v2, *out;
  2445. -ZMAT    *A;
  2446. -complex    alpha;
  2447. -{
  2448. -    /* register    int    j; */
  2449. -    int    i, m, n;
  2450. -    complex    tmp, *v2_ve, *out_ve;
  2451. -    
  2452. -    if ( ! v1 || ! v2 || ! A )
  2453. -    error(E_NULL,"zmv_mltadd");
  2454. -    if ( out == v2 )
  2455. -    error(E_INSITU,"zmv_mltadd");
  2456. -    if ( v1->dim != A->m || v2->dim != A-> n )
  2457. -    error(E_SIZES,"zmv_mltadd");
  2458. -    
  2459. -    tracecatch(out = zv_copy(v1,out),"zmv_mltadd");
  2460. -    
  2461. -    v2_ve = v2->ve;    out_ve = out->ve;
  2462. -    m = A->m;    n = A->n;
  2463. -    
  2464. -    if ( alpha.re == 0.0 && alpha.im == 0.0 )
  2465. -    return out;
  2466. -    
  2467. -    for ( i = 0; i < m; i++ )
  2468. -    {
  2469. -    tmp = __zip__(A->me[i],v2_ve,(int)n,Z_NOCONJ);
  2470. -    out_ve[i].re += alpha.re*tmp.re - alpha.im*tmp.im;
  2471. -    out_ve[i].im += alpha.re*tmp.im + alpha.im*tmp.re;
  2472. -    /**************************************************
  2473. -      A_e = A->me[i];
  2474. -      sum = 0.0;
  2475. -      for ( j = 0; j < n; j++ )
  2476. -      sum += A_e[j]*v2_ve[j];
  2477. -      out_ve[i] = v1->ve[i] + alpha*sum;
  2478. -      **************************************************/
  2479. -    }
  2480. -    
  2481. -    return out;
  2482. -}
  2483. -
  2484. -/* zvm_mltadd -- vector-matrix multiply and add a la zvm_mlt()
  2485. -   -- may not be in situ
  2486. -   -- returns out == v1 + v2*.A */
  2487. -ZVEC    *zvm_mltadd(v1,v2,A,alpha,out)
  2488. -ZVEC    *v1, *v2, *out;
  2489. -ZMAT    *A;
  2490. -complex    alpha;
  2491. -{
  2492. -    int    /* i, */ j, m, n;
  2493. -    complex    tmp, /* *A_e, */ *out_ve;
  2494. -    
  2495. -    if ( ! v1 || ! v2 || ! A )
  2496. -    error(E_NULL,"zvm_mltadd");
  2497. -    if ( v2 == out )
  2498. -    error(E_INSITU,"zvm_mltadd");
  2499. -    if ( v1->dim != A->n || A->m != v2->dim )
  2500. -    error(E_SIZES,"zvm_mltadd");
  2501. -    
  2502. -    tracecatch(out = zv_copy(v1,out),"zvm_mltadd");
  2503. -    
  2504. -    out_ve = out->ve;    m = A->m;    n = A->n;
  2505. -    for ( j = 0; j < m; j++ )
  2506. -    {
  2507. -    /* tmp = zmlt(v2->ve[j],alpha); */
  2508. -    tmp.re =   v2->ve[j].re*alpha.re - v2->ve[j].im*alpha.im;
  2509. -    tmp.im =   v2->ve[j].re*alpha.im + v2->ve[j].im*alpha.re;
  2510. -    if ( tmp.re != 0.0 || tmp.im != 0.0 )
  2511. -        __zmltadd__(out_ve,A->me[j],tmp,(int)n,Z_CONJ);
  2512. -    /**************************************************
  2513. -      A_e = A->me[j];
  2514. -      for ( i = 0; i < n; i++ )
  2515. -      out_ve[i] += A_e[i]*tmp;
  2516. -    **************************************************/
  2517. -    }
  2518. -    
  2519. -    return out;
  2520. -}
  2521. -
  2522. -/* zget_col -- gets a specified column of a matrix; returned as a vector */
  2523. -ZVEC    *zget_col(mat,col,vec)
  2524. -int    col;
  2525. -ZMAT    *mat;
  2526. -ZVEC    *vec;
  2527. -{
  2528. -    u_int    i;
  2529. -
  2530. -    if ( mat==ZMNULL )
  2531. -        error(E_NULL,"zget_col");
  2532. -    if ( col < 0 || col >= mat->n )
  2533. -        error(E_RANGE,"zget_col");
  2534. -    if ( vec==ZVNULL || vec->dim<mat->m )
  2535. -        vec = zv_resize(vec,mat->m);
  2536. -
  2537. -    for ( i=0; i<mat->m; i++ )
  2538. -        vec->ve[i] = mat->me[i][col];
  2539. -
  2540. -    return (vec);
  2541. -}
  2542. -
  2543. -/* zget_row -- gets a specified row of a matrix and retruns it as a vector */
  2544. -ZVEC    *zget_row(mat,row,vec)
  2545. -int    row;
  2546. -ZMAT    *mat;
  2547. -ZVEC    *vec;
  2548. -{
  2549. -    int    /* i, */ lim;
  2550. -
  2551. -    if ( mat==ZMNULL )
  2552. -        error(E_NULL,"zget_row");
  2553. -    if ( row < 0 || row >= mat->m )
  2554. -        error(E_RANGE,"zget_row");
  2555. -    if ( vec==ZVNULL || vec->dim<mat->n )
  2556. -        vec = zv_resize(vec,mat->n);
  2557. -
  2558. -    lim = min(mat->n,vec->dim);
  2559. -
  2560. -    /* for ( i=0; i<mat->n; i++ ) */
  2561. -    /*     vec->ve[i] = mat->me[row][i]; */
  2562. -    MEMCOPY(mat->me[row],vec->ve,lim,complex);
  2563. -
  2564. -    return (vec);
  2565. -}
  2566. -
  2567. -/* zset_col -- sets column of matrix to values given in vec (in situ) */
  2568. -ZMAT    *zset_col(mat,col,vec)
  2569. -ZMAT    *mat;
  2570. -ZVEC    *vec;
  2571. -int    col;
  2572. -{
  2573. -    u_int    i,lim;
  2574. -
  2575. -    if ( mat==ZMNULL || vec==ZVNULL )
  2576. -        error(E_NULL,"zset_col");
  2577. -    if ( col < 0 || col >= mat->n )
  2578. -        error(E_RANGE,"zset_col");
  2579. -    lim = min(mat->m,vec->dim);
  2580. -    for ( i=0; i<lim; i++ )
  2581. -        mat->me[i][col] = vec->ve[i];
  2582. -
  2583. -    return (mat);
  2584. -}
  2585. -
  2586. -/* zset_row -- sets row of matrix to values given in vec (in situ) */
  2587. -ZMAT    *zset_row(mat,row,vec)
  2588. -ZMAT    *mat;
  2589. -ZVEC    *vec;
  2590. -int    row;
  2591. -{
  2592. -    u_int    /* j, */ lim;
  2593. -
  2594. -    if ( mat==ZMNULL || vec==ZVNULL )
  2595. -        error(E_NULL,"zset_row");
  2596. -    if ( row < 0 || row >= mat->m )
  2597. -        error(E_RANGE,"zset_row");
  2598. -    lim = min(mat->n,vec->dim);
  2599. -    /* for ( j=j0; j<lim; j++ ) */
  2600. -    /*     mat->me[row][j] = vec->ve[j]; */
  2601. -    MEMCOPY(vec->ve,mat->me[row],lim,complex);
  2602. -
  2603. -    return (mat);
  2604. -}
  2605. -
  2606. -/* zm_rand -- randomise a complex matrix; uniform in [0,1)+[0,1)*i */
  2607. -ZMAT    *zm_rand(A)
  2608. -ZMAT    *A;
  2609. -{
  2610. -    int        i;
  2611. -
  2612. -    if ( ! A )
  2613. -    error(E_NULL,"zm_rand");
  2614. -
  2615. -    for ( i = 0; i < A->m; i++ )
  2616. -    mrandlist((Real *)(A->me[i]),2*A->n);
  2617. -
  2618. -    return A;
  2619. -}
  2620. //GO.SYSIN DD zmatop.c
  2621. echo znorm.c 1>&2
  2622. sed >znorm.c <<'//GO.SYSIN DD znorm.c' 's/^-//'
  2623. -
  2624. -/**************************************************************************
  2625. -**
  2626. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2627. -**
  2628. -**                 Meschach Library
  2629. -** 
  2630. -** This Meschach Library is provided "as is" without any express 
  2631. -** or implied warranty of any kind with respect to this software. 
  2632. -** In particular the authors shall not be liable for any direct, 
  2633. -** indirect, special, incidental or consequential damages arising 
  2634. -** in any way from use of the software.
  2635. -** 
  2636. -** Everyone is granted permission to copy, modify and redistribute this
  2637. -** Meschach Library, provided:
  2638. -**  1.  All copies contain this copyright notice.
  2639. -**  2.  All modified copies shall carry a notice stating who
  2640. -**      made the last modification and the date of such modification.
  2641. -**  3.  No charge is made for this software or works derived from it.  
  2642. -**      This clause shall not be construed as constraining other software
  2643. -**      distributed on the same medium as this software, nor is a
  2644. -**      distribution fee considered a charge.
  2645. -**
  2646. -***************************************************************************/
  2647. -
  2648. -
  2649. -/*
  2650. -    A collection of functions for computing norms: scaled and unscaled
  2651. -    Complex version
  2652. -*/
  2653. -static    char    rcsid[] = "$Id: znorm.c,v 1.1 1994/01/13 04:21:31 des Exp $";
  2654. -
  2655. -#include    <stdio.h>
  2656. -#include    <math.h>
  2657. -#include    "zmatrix.h"
  2658. -
  2659. -
  2660. -
  2661. -/* _zv_norm1 -- computes (scaled) 1-norms of vectors */
  2662. -double    _zv_norm1(x,scale)
  2663. -ZVEC    *x;
  2664. -VEC    *scale;
  2665. -{
  2666. -    int    i, dim;
  2667. -    Real    s, sum;
  2668. -    
  2669. -    if ( x == ZVNULL )
  2670. -    error(E_NULL,"_zv_norm1");
  2671. -    dim = x->dim;
  2672. -    
  2673. -    sum = 0.0;
  2674. -    if ( scale == VNULL )
  2675. -    for ( i = 0; i < dim; i++ )
  2676. -        sum += zabs(x->ve[i]);
  2677. -    else if ( scale->dim < dim )
  2678. -    error(E_SIZES,"_zv_norm1");
  2679. -    else
  2680. -    for ( i = 0; i < dim; i++ )
  2681. -    {
  2682. -        s = scale->ve[i];
  2683. -        sum += ( s== 0.0 ) ? zabs(x->ve[i]) : zabs(x->ve[i])/fabs(s);
  2684. -    }
  2685. -    
  2686. -    return sum;
  2687. -}
  2688. -
  2689. -/* square -- returns x^2 */
  2690. -/******************************
  2691. -double    square(x)
  2692. -double    x;
  2693. -{    return x*x;    }
  2694. -******************************/
  2695. -
  2696. -#define    square(x)    ((x)*(x))
  2697. -
  2698. -/* _zv_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */
  2699. -double    _zv_norm2(x,scale)
  2700. -ZVEC    *x;
  2701. -VEC    *scale;
  2702. -{
  2703. -    int    i, dim;
  2704. -    Real    s, sum;
  2705. -    
  2706. -    if ( x == ZVNULL )
  2707. -    error(E_NULL,"_zv_norm2");
  2708. -    dim = x->dim;
  2709. -    
  2710. -    sum = 0.0;
  2711. -    if ( scale == VNULL )
  2712. -    for ( i = 0; i < dim; i++ )
  2713. -        sum += square(x->ve[i].re) + square(x->ve[i].im);
  2714. -    else if ( scale->dim < dim )
  2715. -    error(E_SIZES,"_v_norm2");
  2716. -    else
  2717. -    for ( i = 0; i < dim; i++ )
  2718. -    {
  2719. -        s = scale->ve[i];
  2720. -        sum += ( s== 0.0 ) ? square(x->ve[i].re) + square(x->ve[i].im) :
  2721. -        (square(x->ve[i].re) + square(x->ve[i].im))/square(s);
  2722. -    }
  2723. -    
  2724. -    return sqrt(sum);
  2725. -}
  2726. -
  2727. -#define    max(a,b)    ((a) > (b) ? (a) : (b))
  2728. -
  2729. -/* _zv_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */
  2730. -double    _zv_norm_inf(x,scale)
  2731. -ZVEC    *x;
  2732. -VEC    *scale;
  2733. -{
  2734. -    int    i, dim;
  2735. -    Real    s, maxval, tmp;
  2736. -    
  2737. -    if ( x == ZVNULL )
  2738. -    error(E_NULL,"_zv_norm_inf");
  2739. -    dim = x->dim;
  2740. -    
  2741. -    maxval = 0.0;
  2742. -    if ( scale == VNULL )
  2743. -    for ( i = 0; i < dim; i++ )
  2744. -    {
  2745. -        tmp = zabs(x->ve[i]);
  2746. -        maxval = max(maxval,tmp);
  2747. -    }
  2748. -    else if ( scale->dim < dim )
  2749. -    error(E_SIZES,"_zv_norm_inf");
  2750. -    else
  2751. -    for ( i = 0; i < dim; i++ )
  2752. -    {
  2753. -        s = scale->ve[i];
  2754. -        tmp = ( s == 0.0 ) ? zabs(x->ve[i]) : zabs(x->ve[i])/fabs(s);
  2755. -        maxval = max(maxval,tmp);
  2756. -    }
  2757. -    
  2758. -    return maxval;
  2759. -}
  2760. -
  2761. -/* zm_norm1 -- compute matrix 1-norm -- unscaled
  2762. -    -- complex version */
  2763. -double    zm_norm1(A)
  2764. -ZMAT    *A;
  2765. -{
  2766. -    int    i, j, m, n;
  2767. -    Real    maxval, sum;
  2768. -    
  2769. -    if ( A == ZMNULL )
  2770. -    error(E_NULL,"zm_norm1");
  2771. -
  2772. -    m = A->m;    n = A->n;
  2773. -    maxval = 0.0;
  2774. -    
  2775. -    for ( j = 0; j < n; j++ )
  2776. -    {
  2777. -    sum = 0.0;
  2778. -    for ( i = 0; i < m; i ++ )
  2779. -        sum += zabs(A->me[i][j]);
  2780. -    maxval = max(maxval,sum);
  2781. -    }
  2782. -    
  2783. -    return maxval;
  2784. -}
  2785. -
  2786. -/* zm_norm_inf -- compute matrix infinity-norm -- unscaled
  2787. -    -- complex version */
  2788. -double    zm_norm_inf(A)
  2789. -ZMAT    *A;
  2790. -{
  2791. -    int    i, j, m, n;
  2792. -    Real    maxval, sum;
  2793. -    
  2794. -    if ( A == ZMNULL )
  2795. -    error(E_NULL,"zm_norm_inf");
  2796. -    
  2797. -    m = A->m;    n = A->n;
  2798. -    maxval = 0.0;
  2799. -    
  2800. -    for ( i = 0; i < m; i++ )
  2801. -    {
  2802. -    sum = 0.0;
  2803. -    for ( j = 0; j < n; j ++ )
  2804. -        sum += zabs(A->me[i][j]);
  2805. -    maxval = max(maxval,sum);
  2806. -    }
  2807. -    
  2808. -    return maxval;
  2809. -}
  2810. -
  2811. -/* zm_norm_frob -- compute matrix frobenius-norm -- unscaled */
  2812. -double    zm_norm_frob(A)
  2813. -ZMAT    *A;
  2814. -{
  2815. -    int    i, j, m, n;
  2816. -    Real    sum;
  2817. -    
  2818. -    if ( A == ZMNULL )
  2819. -    error(E_NULL,"zm_norm_frob");
  2820. -    
  2821. -    m = A->m;    n = A->n;
  2822. -    sum = 0.0;
  2823. -    
  2824. -    for ( i = 0; i < m; i++ )
  2825. -    for ( j = 0; j < n; j ++ )
  2826. -        sum += square(A->me[i][j].re) + square(A->me[i][j].im);
  2827. -    
  2828. -    return sqrt(sum);
  2829. -}
  2830. -
  2831. //GO.SYSIN DD znorm.c
  2832. echo zfunc.c 1>&2
  2833. sed >zfunc.c <<'//GO.SYSIN DD zfunc.c' 's/^-//'
  2834. -
  2835. -/**************************************************************************
  2836. -**
  2837. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  2838. -**
  2839. -**                 Meschach Library
  2840. -** 
  2841. -** This Meschach Library is provided "as is" without any express 
  2842. -** or implied warranty of any kind with respect to this software. 
  2843. -** In particular the authors shall not be liable for any direct, 
  2844. -** indirect, special, incidental or consequential damages arising 
  2845. -** in any way from use of the software.
  2846. -** 
  2847. -** Everyone is granted permission to copy, modify and redistribute this
  2848. -** Meschach Library, provided:
  2849. -**  1.  All copies contain this copyright notice.
  2850. -**  2.  All modified copies shall carry a notice stating who
  2851. -**      made the last modification and the date of such modification.
  2852. -**  3.  No charge is made for this software or works derived from it.  
  2853. -**      This clause shall not be construed as constraining other software
  2854. -**      distributed on the same medium as this software, nor is a
  2855. -**      distribution fee considered a charge.
  2856. -**
  2857. -***************************************************************************/
  2858. -
  2859. -/*
  2860. -    Elementary functions for complex numbers
  2861. -    -- if not already defined
  2862. -*/
  2863. -
  2864. -#include    <math.h>
  2865. -#include    "zmatrix.h"
  2866. -
  2867. -static char rcsid[] = "$Id: zfunc.c,v 1.1 1994/01/13 04:28:29 des Exp $";
  2868. -
  2869. -#ifndef COMPLEX_H
  2870. -
  2871. -#ifndef zmake
  2872. -/* zmake -- create complex number real + i*imag */
  2873. -complex    zmake(real,imag)
  2874. -double    real, imag;
  2875. -{
  2876. -    complex    w;    /* == real + i*imag */
  2877. -
  2878. -    w.re = real;    w.im = imag;
  2879. -    return w;
  2880. -}
  2881. -#endif
  2882. -
  2883. -#ifndef zneg
  2884. -/* zneg -- returns negative of z */
  2885. -complex    zneg(z)
  2886. -complex    z;
  2887. -{
  2888. -    z.re = - z.re;
  2889. -    z.im = - z.im;
  2890. -
  2891. -    return z;
  2892. -}
  2893. -#endif
  2894. -
  2895. -#ifndef zabs
  2896. -/* zabs -- returns |z| */
  2897. -double    zabs(z)
  2898. -complex    z;
  2899. -{
  2900. -    Real    x, y, tmp;
  2901. -    int        x_expt, y_expt;
  2902. -
  2903. -    /* Note: we must ensure that overflow does not occur! */
  2904. -    x = ( z.re >= 0.0 ) ? z.re : -z.re;
  2905. -    y = ( z.im >= 0.0 ) ? z.im : -z.im;
  2906. -    if ( x < y )
  2907. -    {
  2908. -    tmp = x;
  2909. -    x = y;
  2910. -    y = tmp;
  2911. -    }
  2912. -    x = frexp(x,&x_expt);
  2913. -    y = frexp(y,&y_expt);
  2914. -    y = ldexp(y,y_expt-x_expt);
  2915. -    tmp = sqrt(x*x+y*y);
  2916. -
  2917. -    return ldexp(tmp,x_expt);
  2918. -}
  2919. -#endif
  2920. -
  2921. -#ifndef zadd
  2922. -/* zadd -- returns z1+z2 */
  2923. -complex zadd(z1,z2)
  2924. -complex    z1, z2;
  2925. -{
  2926. -    complex z;
  2927. -
  2928. -    z.re = z1.re + z2.re;
  2929. -    z.im = z1.im + z2.im;
  2930. -
  2931. -    return z;
  2932. -}
  2933. -#endif
  2934. -
  2935. -#ifndef zsub
  2936. -/* zsub -- returns z1-z2 */
  2937. -complex zsub(z1,z2)
  2938. -complex    z1, z2;
  2939. -{
  2940. -    complex z;
  2941. -
  2942. -    z.re = z1.re - z2.re;
  2943. -    z.im = z1.im - z2.im;
  2944. -
  2945. -    return z;
  2946. -}
  2947. -#endif
  2948. -
  2949. -#ifndef zmlt
  2950. -/* zmlt -- returns z1*z2 */
  2951. -complex    zmlt(z1,z2)
  2952. -complex    z1, z2;
  2953. -{
  2954. -    complex z;
  2955. -
  2956. -    z.re = z1.re * z2.re - z1.im * z2.im;
  2957. -    z.im = z1.re * z2.im + z1.im * z2.re;
  2958. -
  2959. -    return z;
  2960. -}
  2961. -#endif
  2962. -
  2963. -#ifndef zinv
  2964. -/* zmlt -- returns 1/z */
  2965. -complex    zinv(z)
  2966. -complex    z;
  2967. -{
  2968. -    Real    x, y, tmp;
  2969. -    int        x_expt, y_expt;
  2970. -
  2971. -    if ( z.re == 0.0 && z.im == 0.0 )
  2972. -    error(E_SING,"zinv");
  2973. -    /* Note: we must ensure that overflow does not occur! */
  2974. -    x = ( z.re >= 0.0 ) ? z.re : -z.re;
  2975. -    y = ( z.im >= 0.0 ) ? z.im : -z.im;
  2976. -    if ( x < y )
  2977. -    {
  2978. -    tmp = x;
  2979. -    x = y;
  2980. -    y = tmp;
  2981. -    }
  2982. -    x = frexp(x,&x_expt);
  2983. -    y = frexp(y,&y_expt);
  2984. -    y = ldexp(y,y_expt-x_expt);
  2985. -
  2986. -    tmp = 1.0/(x*x + y*y);
  2987. -    z.re =  z.re*tmp*ldexp(1.0,-2*x_expt);
  2988. -    z.im = -z.im*tmp*ldexp(1.0,-2*x_expt);
  2989. -
  2990. -    return z;
  2991. -}
  2992. -#endif
  2993. -
  2994. -#ifndef zdiv
  2995. -/* zdiv -- returns z1/z2 */
  2996. -complex    zdiv(z1,z2)
  2997. -complex    z1, z2;
  2998. -{
  2999. -    return zmlt(z1,zinv(z2));
  3000. -}
  3001. -#endif
  3002. -
  3003. -#ifndef zsqrt
  3004. -/* zsqrt -- returns sqrt(z); uses branch with Re sqrt(z) >= 0 */
  3005. -complex    zsqrt(z)
  3006. -complex    z;
  3007. -{
  3008. -    complex    w;    /* == sqrt(z) at end */
  3009. -    Real    alpha;
  3010. -
  3011. -    alpha = sqrt(0.5*(fabs(z.re) + zabs(z)));
  3012. -    if ( z.re >= 0.0 )
  3013. -    {
  3014. -    w.re = alpha;
  3015. -    w.im = z.im / (2.0*alpha);
  3016. -    }
  3017. -    else
  3018. -    {
  3019. -    w.re = fabs(z.im)/(2.0*alpha);
  3020. -    w.im = ( z.im >= 0 ) ? alpha : - alpha;
  3021. -    }
  3022. -
  3023. -    return w;
  3024. -}
  3025. -#endif
  3026. -
  3027. -#ifndef    zexp
  3028. -/* zexp -- returns exp(z) */
  3029. -complex    zexp(z)
  3030. -complex    z;
  3031. -{
  3032. -    complex    w;    /* == exp(z) at end */
  3033. -    Real    r;
  3034. -
  3035. -    r = exp(z.re);
  3036. -    w.re = r*cos(z.im);
  3037. -    w.im = r*sin(z.im);
  3038. -
  3039. -    return w;
  3040. -}
  3041. -#endif
  3042. -
  3043. -#ifndef    zlog
  3044. -/* zlog -- returns log(z); uses principal branch with -pi <= Im log(z) <= pi */
  3045. -complex    zlog(z)
  3046. -complex    z;
  3047. -{
  3048. -    complex    w;    /* == log(z) at end */
  3049. -
  3050. -    w.re = log(zabs(z));
  3051. -    w.im = atan2(z.im,z.re);
  3052. -
  3053. -    return w;
  3054. -}
  3055. -#endif
  3056. -
  3057. -#ifndef zconj
  3058. -complex    zconj(z)
  3059. -complex    z;
  3060. -{
  3061. -    complex    w;    /* == conj(z) */
  3062. -
  3063. -    w.re =   z.re;
  3064. -    w.im = - z.im;
  3065. -    return w;
  3066. -}
  3067. -#endif
  3068. -
  3069. -#endif
  3070. //GO.SYSIN DD zfunc.c
  3071. echo zlufctr.c 1>&2
  3072. sed >zlufctr.c <<'//GO.SYSIN DD zlufctr.c' 's/^-//'
  3073. -
  3074. -/**************************************************************************
  3075. -**
  3076. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3077. -**
  3078. -**                 Meschach Library
  3079. -** 
  3080. -** This Meschach Library is provided "as is" without any express 
  3081. -** or implied warranty of any kind with respect to this software. 
  3082. -** In particular the authors shall not be liable for any direct, 
  3083. -** indirect, special, incidental or consequential damages arising 
  3084. -** in any way from use of the software.
  3085. -** 
  3086. -** Everyone is granted permission to copy, modify and redistribute this
  3087. -** Meschach Library, provided:
  3088. -**  1.  All copies contain this copyright notice.
  3089. -**  2.  All modified copies shall carry a notice stating who
  3090. -**      made the last modification and the date of such modification.
  3091. -**  3.  No charge is made for this software or works derived from it.  
  3092. -**      This clause shall not be construed as constraining other software
  3093. -**      distributed on the same medium as this software, nor is a
  3094. -**      distribution fee considered a charge.
  3095. -**
  3096. -***************************************************************************/
  3097. -
  3098. -/*
  3099. -    Matrix factorisation routines to work with the other matrix files.
  3100. -    Complex version
  3101. -*/
  3102. -
  3103. -static    char    rcsid[] = "$Id: zlufctr.c,v 1.1 1994/01/13 04:26:20 des Exp $";
  3104. -
  3105. -#include    <stdio.h>
  3106. -#include    <math.h>
  3107. -#include    "zmatrix.h"
  3108. -#include        "zmatrix2.h"
  3109. -
  3110. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  3111. -
  3112. -
  3113. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  3114. -
  3115. -/* zLUfactor -- Gaussian elimination with scaled partial pivoting
  3116. -        -- Note: returns LU matrix which is A */
  3117. -ZMAT    *zLUfactor(A,pivot)
  3118. -ZMAT    *A;
  3119. -PERM    *pivot;
  3120. -{
  3121. -    u_int    i, j, k, k_max, m, n;
  3122. -    int    i_max;
  3123. -    Real    dtemp, max1;
  3124. -    complex    **A_v, *A_piv, *A_row, temp;
  3125. -    static    VEC    *scale = VNULL;
  3126. -
  3127. -    if ( A==ZMNULL || pivot==PNULL )
  3128. -        error(E_NULL,"zLUfactor");
  3129. -    if ( pivot->size != A->m )
  3130. -        error(E_SIZES,"zLUfactor");
  3131. -    m = A->m;    n = A->n;
  3132. -    scale = v_resize(scale,A->m);
  3133. -    MEM_STAT_REG(scale,TYPE_VEC);
  3134. -    A_v = A->me;
  3135. -
  3136. -    /* initialise pivot with identity permutation */
  3137. -    for ( i=0; i<m; i++ )
  3138. -        pivot->pe[i] = i;
  3139. -
  3140. -    /* set scale parameters */
  3141. -    for ( i=0; i<m; i++ )
  3142. -    {
  3143. -        max1 = 0.0;
  3144. -        for ( j=0; j<n; j++ )
  3145. -        {
  3146. -            dtemp = zabs(A_v[i][j]);
  3147. -            max1 = max(max1,dtemp);
  3148. -        }
  3149. -        scale->ve[i] = max1;
  3150. -    }
  3151. -
  3152. -    /* main loop */
  3153. -    k_max = min(m,n)-1;
  3154. -    for ( k=0; k<k_max; k++ )
  3155. -    {
  3156. -        /* find best pivot row */
  3157. -        max1 = 0.0;    i_max = -1;
  3158. -        for ( i=k; i<m; i++ )
  3159. -        if ( scale->ve[i] > 0.0 )
  3160. -        {
  3161. -            dtemp = zabs(A_v[i][k])/scale->ve[i];
  3162. -            if ( dtemp > max1 )
  3163. -            { max1 = dtemp;    i_max = i;    }
  3164. -        }
  3165. -        
  3166. -        /* if no pivot then ignore column k... */
  3167. -        if ( i_max == -1 )
  3168. -        continue;
  3169. -
  3170. -        /* do we pivot ? */
  3171. -        if ( i_max != k )    /* yes we do... */
  3172. -        {
  3173. -        px_transp(pivot,i_max,k);
  3174. -        for ( j=0; j<n; j++ )
  3175. -        {
  3176. -            temp = A_v[i_max][j];
  3177. -            A_v[i_max][j] = A_v[k][j];
  3178. -            A_v[k][j] = temp;
  3179. -        }
  3180. -        }
  3181. -        
  3182. -        /* row operations */
  3183. -        for ( i=k+1; i<m; i++ )    /* for each row do... */
  3184. -        {    /* Note: divide by zero should never happen */
  3185. -        temp = A_v[i][k] = zdiv(A_v[i][k],A_v[k][k]);
  3186. -        A_piv = &(A_v[k][k+1]);
  3187. -        A_row = &(A_v[i][k+1]);
  3188. -        temp.re = - temp.re;
  3189. -        temp.im = - temp.im;
  3190. -        if ( k+1 < n )
  3191. -            __zmltadd__(A_row,A_piv,temp,(int)(n-(k+1)),Z_NOCONJ);
  3192. -        /*********************************************
  3193. -          for ( j=k+1; j<n; j++ )
  3194. -          A_v[i][j] -= temp*A_v[k][j];
  3195. -          (*A_row++) -= temp*(*A_piv++);
  3196. -        *********************************************/
  3197. -        }
  3198. -    }
  3199. -
  3200. -    return A;
  3201. -}
  3202. -
  3203. -
  3204. -/* zLUsolve -- given an LU factorisation in A, solve Ax=b */
  3205. -ZVEC    *zLUsolve(A,pivot,b,x)
  3206. -ZMAT    *A;
  3207. -PERM    *pivot;
  3208. -ZVEC    *b,*x;
  3209. -{
  3210. -    if ( A==ZMNULL || b==ZVNULL || pivot==PNULL )
  3211. -        error(E_NULL,"zLUsolve");
  3212. -    if ( A->m != A->n || A->n != b->dim )
  3213. -        error(E_SIZES,"zLUsolve");
  3214. -
  3215. -    x = px_zvec(pivot,b,x);    /* x := P.b */
  3216. -    zLsolve(A,x,x,1.0);    /* implicit diagonal = 1 */
  3217. -    zUsolve(A,x,x,0.0);    /* explicit diagonal */
  3218. -
  3219. -    return (x);
  3220. -}
  3221. -
  3222. -/* zLUAsolve -- given an LU factorisation in A, solve A^*.x=b */
  3223. -ZVEC    *zLUAsolve(LU,pivot,b,x)
  3224. -ZMAT    *LU;
  3225. -PERM    *pivot;
  3226. -ZVEC    *b,*x;
  3227. -{
  3228. -    if ( ! LU || ! b || ! pivot )
  3229. -        error(E_NULL,"zLUAsolve");
  3230. -    if ( LU->m != LU->n || LU->n != b->dim )
  3231. -        error(E_SIZES,"zLUAsolve");
  3232. -
  3233. -    x = zv_copy(b,x);
  3234. -    zUAsolve(LU,x,x,0.0);    /* explicit diagonal */
  3235. -    zLAsolve(LU,x,x,1.0);    /* implicit diagonal = 1 */
  3236. -    pxinv_zvec(pivot,x,x);    /* x := P^*.x */
  3237. -
  3238. -    return (x);
  3239. -}
  3240. -
  3241. -/* zm_inverse -- returns inverse of A, provided A is not too rank deficient
  3242. -    -- uses LU factorisation */
  3243. -ZMAT    *zm_inverse(A,out)
  3244. -ZMAT    *A, *out;
  3245. -{
  3246. -    int    i;
  3247. -    ZVEC    *tmp, *tmp2;
  3248. -    ZMAT    *A_cp;
  3249. -    PERM    *pivot;
  3250. -
  3251. -    if ( ! A )
  3252. -        error(E_NULL,"zm_inverse");
  3253. -    if ( A->m != A->n )
  3254. -        error(E_SQUARE,"zm_inverse");
  3255. -    if ( ! out || out->m < A->m || out->n < A->n )
  3256. -        out = zm_resize(out,A->m,A->n);
  3257. -
  3258. -    A_cp = zm_copy(A,ZMNULL);
  3259. -    tmp = zv_get(A->m);
  3260. -    tmp2 = zv_get(A->m);
  3261. -    pivot = px_get(A->m);
  3262. -    tracecatch(zLUfactor(A_cp,pivot),"zm_inverse");
  3263. -    for ( i = 0; i < A->n; i++ )
  3264. -    {
  3265. -        zv_zero(tmp);
  3266. -        tmp->ve[i].re = 1.0;
  3267. -        tmp->ve[i].im = 0.0;
  3268. -        tracecatch(zLUsolve(A_cp,pivot,tmp,tmp2),"m_inverse");
  3269. -        zset_col(out,i,tmp2);
  3270. -    }
  3271. -
  3272. -    ZM_FREE(A_cp);
  3273. -    ZV_FREE(tmp);    ZV_FREE(tmp2);
  3274. -    PX_FREE(pivot);
  3275. -
  3276. -    return out;
  3277. -}
  3278. -
  3279. -/* zLUcondest -- returns an estimate of the condition number of LU given the
  3280. -    LU factorisation in compact form */
  3281. -double    zLUcondest(LU,pivot)
  3282. -ZMAT    *LU;
  3283. -PERM    *pivot;
  3284. -{
  3285. -    static    ZVEC    *y = ZVNULL, *z = ZVNULL;
  3286. -    Real    cond_est, L_norm, U_norm, norm, sn_inv;
  3287. -    complex    sum;
  3288. -    int        i, j, n;
  3289. -
  3290. -    if ( ! LU || ! pivot )
  3291. -    error(E_NULL,"zLUcondest");
  3292. -    if ( LU->m != LU->n )
  3293. -    error(E_SQUARE,"zLUcondest");
  3294. -    if ( LU->n != pivot->size )
  3295. -    error(E_SIZES,"zLUcondest");
  3296. -
  3297. -    n = LU->n;
  3298. -    y = zv_resize(y,n);
  3299. -    z = zv_resize(z,n);
  3300. -    MEM_STAT_REG(y,TYPE_ZVEC);
  3301. -    MEM_STAT_REG(z,TYPE_ZVEC);
  3302. -
  3303. -    cond_est = 0.0;        /* should never be returned */
  3304. -
  3305. -    for ( i = 0; i < n; i++ )
  3306. -    {
  3307. -    sum.re = 1.0;
  3308. -    sum.im = 0.0;
  3309. -    for ( j = 0; j < i; j++ )
  3310. -        /* sum -= LU->me[j][i]*y->ve[j]; */
  3311. -        sum = zsub(sum,zmlt(LU->me[j][i],y->ve[j]));
  3312. -    /* sum -= (sum < 0.0) ? 1.0 : -1.0; */
  3313. -    sn_inv = 1.0 / zabs(sum);
  3314. -    sum.re += sum.re * sn_inv;
  3315. -    sum.im += sum.im * sn_inv;
  3316. -    if ( is_zero(LU->me[i][i]) )
  3317. -        return HUGE;
  3318. -    /* y->ve[i] = sum / LU->me[i][i]; */
  3319. -    y->ve[i] = zdiv(sum,LU->me[i][i]);
  3320. -    }
  3321. -
  3322. -    zLAsolve(LU,y,y,1.0);
  3323. -    zLUsolve(LU,pivot,y,z);
  3324. -
  3325. -    /* now estimate norm of A (even though it is not directly available) */
  3326. -    /* actually computes ||L||_inf.||U||_inf */
  3327. -    U_norm = 0.0;
  3328. -    for ( i = 0; i < n; i++ )
  3329. -    {
  3330. -    norm = 0.0;
  3331. -    for ( j = i; j < n; j++ )
  3332. -        norm += zabs(LU->me[i][j]);
  3333. -    if ( norm > U_norm )
  3334. -        U_norm = norm;
  3335. -    }
  3336. -    L_norm = 0.0;
  3337. -    for ( i = 0; i < n; i++ )
  3338. -    {
  3339. -    norm = 1.0;
  3340. -    for ( j = 0; j < i; j++ )
  3341. -        norm += zabs(LU->me[i][j]);
  3342. -    if ( norm > L_norm )
  3343. -        L_norm = norm;
  3344. -    }
  3345. -
  3346. -    tracecatch(cond_est = U_norm*L_norm*zv_norm_inf(z)/zv_norm_inf(y),
  3347. -           "LUcondest");
  3348. -
  3349. -    return cond_est;
  3350. -}
  3351. //GO.SYSIN DD zlufctr.c
  3352. echo zsolve.c 1>&2
  3353. sed >zsolve.c <<'//GO.SYSIN DD zsolve.c' 's/^-//'
  3354. -
  3355. -/**************************************************************************
  3356. -**
  3357. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3358. -**
  3359. -**                 Meschach Library
  3360. -** 
  3361. -** This Meschach Library is provided "as is" without any express 
  3362. -** or implied warranty of any kind with respect to this software. 
  3363. -** In particular the authors shall not be liable for any direct, 
  3364. -** indirect, special, incidental or consequential damages arising 
  3365. -** in any way from use of the software.
  3366. -** 
  3367. -** Everyone is granted permission to copy, modify and redistribute this
  3368. -** Meschach Library, provided:
  3369. -**  1.  All copies contain this copyright notice.
  3370. -**  2.  All modified copies shall carry a notice stating who
  3371. -**      made the last modification and the date of such modification.
  3372. -**  3.  No charge is made for this software or works derived from it.  
  3373. -**      This clause shall not be construed as constraining other software
  3374. -**      distributed on the same medium as this software, nor is a
  3375. -**      distribution fee considered a charge.
  3376. -**
  3377. -***************************************************************************/
  3378. -
  3379. -
  3380. -/*
  3381. -    Matrix factorisation routines to work with the other matrix files.
  3382. -    Complex case
  3383. -*/
  3384. -
  3385. -static    char    rcsid[] = "$Id: zsolve.c,v 1.1 1994/01/13 04:20:33 des Exp $";
  3386. -
  3387. -#include    <stdio.h>
  3388. -#include    <math.h>
  3389. -#include        "zmatrix2.h"
  3390. -
  3391. -
  3392. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0 )
  3393. -
  3394. -/* Most matrix factorisation routines are in-situ unless otherwise specified */
  3395. -
  3396. -/* zUsolve -- back substitution with optional over-riding diagonal
  3397. -        -- can be in-situ but doesn't need to be */
  3398. -ZVEC    *zUsolve(matrix,b,out,diag)
  3399. -ZMAT    *matrix;
  3400. -ZVEC    *b, *out;
  3401. -double    diag;
  3402. -{
  3403. -    u_int    dim /* , j */;
  3404. -    int    i, i_lim;
  3405. -    complex    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  3406. -    
  3407. -    if ( matrix==ZMNULL || b==ZVNULL )
  3408. -    error(E_NULL,"zUsolve");
  3409. -    dim = min(matrix->m,matrix->n);
  3410. -    if ( b->dim < dim )
  3411. -    error(E_SIZES,"zUsolve");
  3412. -    if ( out==ZVNULL || out->dim < dim )
  3413. -    out = zv_resize(out,matrix->n);
  3414. -    mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  3415. -    
  3416. -    for ( i=dim-1; i>=0; i-- )
  3417. -    if ( ! is_zero(b_ent[i]) )
  3418. -        break;
  3419. -    else
  3420. -        out_ent[i].re = out_ent[i].im = 0.0;
  3421. -    i_lim = i;
  3422. -    
  3423. -    for ( i = i_lim; i>=0; i-- )
  3424. -    {
  3425. -    sum = b_ent[i];
  3426. -    mat_row = &(mat_ent[i][i+1]);
  3427. -    out_col = &(out_ent[i+1]);
  3428. -    sum = zsub(sum,__zip__(mat_row,out_col,i_lim-i,Z_NOCONJ));
  3429. -    /******************************************************
  3430. -      for ( j=i+1; j<=i_lim; j++ )
  3431. -      sum -= mat_ent[i][j]*out_ent[j];
  3432. -      sum -= (*mat_row++)*(*out_col++);
  3433. -    ******************************************************/
  3434. -    if ( diag == 0.0 )
  3435. -    {
  3436. -        if ( is_zero(mat_ent[i][i]) )
  3437. -        error(E_SING,"zUsolve");
  3438. -        else
  3439. -        /* out_ent[i] = sum/mat_ent[i][i]; */
  3440. -        out_ent[i] = zdiv(sum,mat_ent[i][i]);
  3441. -    }
  3442. -    else
  3443. -    {
  3444. -        /* out_ent[i] = sum/diag; */
  3445. -        out_ent[i].re = sum.re / diag;
  3446. -        out_ent[i].im = sum.im / diag;
  3447. -    }
  3448. -    }
  3449. -    
  3450. -    return (out);
  3451. -}
  3452. -
  3453. -/* zLsolve -- forward elimination with (optional) default diagonal value */
  3454. -ZVEC    *zLsolve(matrix,b,out,diag)
  3455. -ZMAT    *matrix;
  3456. -ZVEC    *b,*out;
  3457. -double    diag;
  3458. -{
  3459. -    u_int    dim, i, i_lim /* , j */;
  3460. -    complex    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  3461. -    
  3462. -    if ( matrix==ZMNULL || b==ZVNULL )
  3463. -    error(E_NULL,"zLsolve");
  3464. -    dim = min(matrix->m,matrix->n);
  3465. -    if ( b->dim < dim )
  3466. -    error(E_SIZES,"zLsolve");
  3467. -    if ( out==ZVNULL || out->dim < dim )
  3468. -    out = zv_resize(out,matrix->n);
  3469. -    mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  3470. -    
  3471. -    for ( i=0; i<dim; i++ )
  3472. -    if ( ! is_zero(b_ent[i]) )
  3473. -        break;
  3474. -    else
  3475. -        out_ent[i].re = out_ent[i].im = 0.0;
  3476. -    i_lim = i;
  3477. -    
  3478. -    for ( i = i_lim; i<dim; i++ )
  3479. -    {
  3480. -    sum = b_ent[i];
  3481. -    mat_row = &(mat_ent[i][i_lim]);
  3482. -    out_col = &(out_ent[i_lim]);
  3483. -    sum = zsub(sum,__zip__(mat_row,out_col,(int)(i-i_lim),Z_NOCONJ));
  3484. -    /*****************************************************
  3485. -      for ( j=i_lim; j<i; j++ )
  3486. -      sum -= mat_ent[i][j]*out_ent[j];
  3487. -      sum -= (*mat_row++)*(*out_col++);
  3488. -    ******************************************************/
  3489. -    if ( diag == 0.0 )
  3490. -    {
  3491. -        if ( is_zero(mat_ent[i][i]) )
  3492. -        error(E_SING,"zLsolve");
  3493. -        else
  3494. -        out_ent[i] = zdiv(sum,mat_ent[i][i]);
  3495. -    }
  3496. -    else
  3497. -    {
  3498. -        out_ent[i].re = sum.re / diag;
  3499. -        out_ent[i].im = sum.im / diag;
  3500. -    }
  3501. -    }
  3502. -    
  3503. -    return (out);
  3504. -}
  3505. -
  3506. -
  3507. -/* zUAsolve -- forward elimination with (optional) default diagonal value
  3508. -        using UPPER triangular part of matrix */
  3509. -ZVEC    *zUAsolve(U,b,out,diag)
  3510. -ZMAT    *U;
  3511. -ZVEC    *b,*out;
  3512. -double    diag;
  3513. -{
  3514. -    u_int    dim, i, i_lim /* , j */;
  3515. -    complex    **U_me, *b_ve, *out_ve, tmp;
  3516. -    Real    invdiag;
  3517. -    
  3518. -    if ( ! U || ! b )
  3519. -    error(E_NULL,"zUAsolve");
  3520. -    dim = min(U->m,U->n);
  3521. -    if ( b->dim < dim )
  3522. -    error(E_SIZES,"zUAsolve");
  3523. -    out = zv_resize(out,U->n);
  3524. -    U_me = U->me;    b_ve = b->ve;    out_ve = out->ve;
  3525. -    
  3526. -    for ( i=0; i<dim; i++ )
  3527. -    if ( ! is_zero(b_ve[i]) )
  3528. -        break;
  3529. -    else
  3530. -        out_ve[i].re = out_ve[i].im = 0.0;
  3531. -    i_lim = i;
  3532. -    if ( b != out )
  3533. -    {
  3534. -    __zzero__(out_ve,out->dim);
  3535. -    /* MEM_COPY(&(b_ve[i_lim]),&(out_ve[i_lim]),
  3536. -       (dim-i_lim)*sizeof(complex)); */
  3537. -    MEMCOPY(&(b_ve[i_lim]),&(out_ve[i_lim]),dim-i_lim,complex);
  3538. -    }
  3539. -
  3540. -    if ( diag == 0.0 )
  3541. -    {
  3542. -    for (    ; i<dim; i++ )
  3543. -    {
  3544. -        tmp = zconj(U_me[i][i]);
  3545. -        if ( is_zero(tmp) )
  3546. -        error(E_SING,"zUAsolve");
  3547. -        /* out_ve[i] /= tmp; */
  3548. -        out_ve[i] = zdiv(out_ve[i],tmp);
  3549. -        tmp.re = - out_ve[i].re;
  3550. -        tmp.im = - out_ve[i].im;
  3551. -        __zmltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),tmp,dim-i-1,Z_CONJ);
  3552. -    }
  3553. -    }
  3554. -    else
  3555. -    {
  3556. -    invdiag = 1.0/diag;
  3557. -    for (    ; i<dim; i++ )
  3558. -    {
  3559. -        out_ve[i].re *= invdiag;
  3560. -        out_ve[i].im *= invdiag;
  3561. -        tmp.re = - out_ve[i].re;
  3562. -        tmp.im = - out_ve[i].im;
  3563. -        __zmltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),tmp,dim-i-1,Z_CONJ);
  3564. -    }
  3565. -    }
  3566. -    return (out);
  3567. -}
  3568. -
  3569. -/* zDsolve -- solves Dx=b where D is the diagonal of A -- may be in-situ */
  3570. -ZVEC    *zDsolve(A,b,x)
  3571. -ZMAT    *A;
  3572. -ZVEC    *b,*x;
  3573. -{
  3574. -    u_int    dim, i;
  3575. -    
  3576. -    if ( ! A || ! b )
  3577. -    error(E_NULL,"zDsolve");
  3578. -    dim = min(A->m,A->n);
  3579. -    if ( b->dim < dim )
  3580. -    error(E_SIZES,"zDsolve");
  3581. -    x = zv_resize(x,A->n);
  3582. -    
  3583. -    dim = b->dim;
  3584. -    for ( i=0; i<dim; i++ )
  3585. -    if ( is_zero(A->me[i][i]) )
  3586. -        error(E_SING,"zDsolve");
  3587. -    else
  3588. -        x->ve[i] = zdiv(b->ve[i],A->me[i][i]);
  3589. -    
  3590. -    return (x);
  3591. -}
  3592. -
  3593. -/* zLAsolve -- back substitution with optional over-riding diagonal
  3594. -        using the LOWER triangular part of matrix
  3595. -        -- can be in-situ but doesn't need to be */
  3596. -ZVEC    *zLAsolve(L,b,out,diag)
  3597. -ZMAT    *L;
  3598. -ZVEC    *b, *out;
  3599. -double    diag;
  3600. -{
  3601. -    u_int    dim;
  3602. -    int        i, i_lim;
  3603. -    complex    **L_me, *b_ve, *out_ve, tmp;
  3604. -    Real    invdiag;
  3605. -    
  3606. -    if ( ! L || ! b )
  3607. -    error(E_NULL,"zLAsolve");
  3608. -    dim = min(L->m,L->n);
  3609. -    if ( b->dim < dim )
  3610. -    error(E_SIZES,"zLAsolve");
  3611. -    out = zv_resize(out,L->n);
  3612. -    L_me = L->me;    b_ve = b->ve;    out_ve = out->ve;
  3613. -    
  3614. -    for ( i=dim-1; i>=0; i-- )
  3615. -    if ( ! is_zero(b_ve[i]) )
  3616. -        break;
  3617. -    i_lim = i;
  3618. -
  3619. -    if ( b != out )
  3620. -    {
  3621. -    __zzero__(out_ve,out->dim);
  3622. -    /* MEM_COPY(b_ve,out_ve,(i_lim+1)*sizeof(complex)); */
  3623. -    MEMCOPY(b_ve,out_ve,i_lim+1,complex);
  3624. -    }
  3625. -
  3626. -    if ( diag == 0.0 )
  3627. -    {
  3628. -    for (        ; i>=0; i-- )
  3629. -    {
  3630. -        tmp = zconj(L_me[i][i]);
  3631. -        if ( is_zero(tmp) )
  3632. -        error(E_SING,"zLAsolve");
  3633. -        out_ve[i] = zdiv(out_ve[i],tmp);
  3634. -        tmp.re = - out_ve[i].re;
  3635. -        tmp.im = - out_ve[i].im;
  3636. -        __zmltadd__(out_ve,L_me[i],tmp,i,Z_CONJ);
  3637. -    }
  3638. -    }
  3639. -    else
  3640. -    {
  3641. -    invdiag = 1.0/diag;
  3642. -    for (        ; i>=0; i-- )
  3643. -    {
  3644. -        out_ve[i].re *= invdiag;
  3645. -        out_ve[i].im *= invdiag;
  3646. -        tmp.re = - out_ve[i].re;
  3647. -        tmp.im = - out_ve[i].im;
  3648. -        __zmltadd__(out_ve,L_me[i],tmp,i,Z_CONJ);
  3649. -    }
  3650. -    }
  3651. -    
  3652. -    return (out);
  3653. -}
  3654. //GO.SYSIN DD zsolve.c
  3655. echo zmatlab.c 1>&2
  3656. sed >zmatlab.c <<'//GO.SYSIN DD zmatlab.c' 's/^-//'
  3657. -
  3658. -/**************************************************************************
  3659. -**
  3660. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3661. -**
  3662. -**                 Meschach Library
  3663. -** 
  3664. -** This Meschach Library is provided "as is" without any express 
  3665. -** or implied warranty of any kind with respect to this software. 
  3666. -** In particular the authors shall not be liable for any direct, 
  3667. -** indirect, special, incidental or consequential damages arising 
  3668. -** in any way from use of the software.
  3669. -** 
  3670. -** Everyone is granted permission to copy, modify and redistribute this
  3671. -** Meschach Library, provided:
  3672. -**  1.  All copies contain this copyright notice.
  3673. -**  2.  All modified copies shall carry a notice stating who
  3674. -**      made the last modification and the date of such modification.
  3675. -**  3.  No charge is made for this software or works derived from it.  
  3676. -**      This clause shall not be construed as constraining other software
  3677. -**      distributed on the same medium as this software, nor is a
  3678. -**      distribution fee considered a charge.
  3679. -**
  3680. -***************************************************************************/
  3681. -
  3682. -
  3683. -/*
  3684. -    This file contains routines for import/exporting complex data
  3685. -    to/from MATLAB. The main routines are:
  3686. -            ZMAT *zm_save(FILE *fp,ZMAT *A,char *name)
  3687. -            ZVEC *zv_save(FILE *fp,ZVEC *x,char *name)
  3688. -            complex z_save(FILE *fp,complex z,char *name)
  3689. -            ZMAT *zm_load(FILE *fp,char **name)
  3690. -*/
  3691. -
  3692. -#include        <stdio.h>
  3693. -#include        "zmatrix.h"
  3694. -#include    "matlab.h"
  3695. -
  3696. -static char rcsid[] = "$Id: zmatlab.c,v 1.1 1994/01/13 04:24:57 des Exp $";
  3697. -
  3698. -/* zm_save -- save matrix in ".mat" file for MATLAB
  3699. -   -- returns matrix to be saved */
  3700. -ZMAT    *zm_save(fp,A,name)
  3701. -FILE    *fp;
  3702. -ZMAT    *A;
  3703. -char    *name;
  3704. -{
  3705. -    int     i, j;
  3706. -    matlab  mat;
  3707. -    
  3708. -    if ( ! A )
  3709. -    error(E_NULL,"zm_save");
  3710. -    
  3711. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  3712. -    mat.m = A->m;
  3713. -    mat.n = A->n;
  3714. -    mat.imag = TRUE;
  3715. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  3716. -    
  3717. -    /* write header */
  3718. -    fwrite(&mat,sizeof(matlab),1,fp);
  3719. -    /* write name */
  3720. -    if ( name == (char *)NULL )
  3721. -    fwrite("",sizeof(char),1,fp);
  3722. -    else
  3723. -    fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  3724. -    /* write actual data */
  3725. -    for ( i = 0; i < A->m; i++ )
  3726. -    for ( j = 0; j < A->n; j++ )
  3727. -        fwrite(&(A->me[i][j].re),sizeof(Real),1,fp);
  3728. -    for ( i = 0; i < A->m; i++ )
  3729. -    for ( j = 0; j < A->n; j++ )
  3730. -        fwrite(&(A->me[i][j].im),sizeof(Real),1,fp);
  3731. -    
  3732. -    return A;
  3733. -}
  3734. -
  3735. -
  3736. -/* zv_save -- save vector in ".mat" file for MATLAB
  3737. -   -- saves it as a row vector
  3738. -   -- returns vector to be saved */
  3739. -ZVEC    *zv_save(fp,x,name)
  3740. -FILE    *fp;
  3741. -ZVEC    *x;
  3742. -char    *name;
  3743. -{
  3744. -    int    i;
  3745. -    matlab  mat;
  3746. -    
  3747. -    if ( ! x )
  3748. -    error(E_NULL,"zv_save");
  3749. -    
  3750. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  3751. -    mat.m = x->dim;
  3752. -    mat.n = 1;
  3753. -    mat.imag = TRUE;
  3754. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  3755. -    
  3756. -    /* write header */
  3757. -    fwrite(&mat,sizeof(matlab),1,fp);
  3758. -    /* write name */
  3759. -    if ( name == (char *)NULL )
  3760. -    fwrite("",sizeof(char),1,fp);
  3761. -    else
  3762. -    fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  3763. -    /* write actual data */
  3764. -    for ( i = 0; i < x->dim; i++ )
  3765. -    fwrite(&(x->ve[i].re),sizeof(Real),1,fp);
  3766. -    for ( i = 0; i < x->dim; i++ )
  3767. -    fwrite(&(x->ve[i].im),sizeof(Real),1,fp);
  3768. -    
  3769. -    return x;
  3770. -}
  3771. -
  3772. -/* z_save -- saves complex number in ".mat" file for MATLAB
  3773. -    -- returns complex number to be saved */
  3774. -complex    z_save(fp,z,name)
  3775. -FILE    *fp;
  3776. -complex    z;
  3777. -char    *name;
  3778. -{
  3779. -    matlab  mat;
  3780. -    
  3781. -    mat.type = 1000*MACH_ID + 100*ORDER + 10*PRECISION + 0;
  3782. -    mat.m = 1;
  3783. -    mat.n = 1;
  3784. -    mat.imag = TRUE;
  3785. -    mat.namlen = (name == (char *)NULL) ? 1 : strlen(name)+1;
  3786. -    
  3787. -    /* write header */
  3788. -    fwrite(&mat,sizeof(matlab),1,fp);
  3789. -    /* write name */
  3790. -    if ( name == (char *)NULL )
  3791. -    fwrite("",sizeof(char),1,fp);
  3792. -    else
  3793. -    fwrite(name,sizeof(char),(int)(mat.namlen),fp);
  3794. -    /* write actual data */
  3795. -    fwrite(&z,sizeof(complex),1,fp);
  3796. -    
  3797. -    return z;
  3798. -}
  3799. -
  3800. -
  3801. -
  3802. -/* zm_load -- loads in a ".mat" file variable as produced by MATLAB
  3803. -   -- matrix returned; imaginary parts ignored */
  3804. -ZMAT    *zm_load(fp,name)
  3805. -FILE    *fp;
  3806. -char    **name;
  3807. -{
  3808. -    ZMAT     *A;
  3809. -    int     i;
  3810. -    int     m_flag, o_flag, p_flag, t_flag;
  3811. -    float   f_temp;
  3812. -    double  d_temp;
  3813. -    matlab  mat;
  3814. -    
  3815. -    if ( fread(&mat,sizeof(matlab),1,fp) != 1 )
  3816. -    error(E_FORMAT,"zm_load");
  3817. -    if ( mat.type >= 10000 )    /* don't load a sparse matrix! */
  3818. -    error(E_FORMAT,"zm_load");
  3819. -    m_flag = (mat.type/1000) % 10;
  3820. -    o_flag = (mat.type/100) % 10;
  3821. -    p_flag = (mat.type/10) % 10;
  3822. -    t_flag = (mat.type) % 10;
  3823. -    if ( m_flag != MACH_ID )
  3824. -    error(E_FORMAT,"zm_load");
  3825. -    if ( t_flag != 0 )
  3826. -    error(E_FORMAT,"zm_load");
  3827. -    if ( p_flag != DOUBLE_PREC && p_flag != SINGLE_PREC )
  3828. -    error(E_FORMAT,"zm_load");
  3829. -    *name = (char *)malloc((unsigned)(mat.namlen)+1);
  3830. -    if ( fread(*name,sizeof(char),(unsigned)(mat.namlen),fp) == 0 )
  3831. -    error(E_FORMAT,"zm_load");
  3832. -    A = zm_get((unsigned)(mat.m),(unsigned)(mat.n));
  3833. -    for ( i = 0; i < A->m*A->n; i++ )
  3834. -    {
  3835. -    if ( p_flag == DOUBLE_PREC )
  3836. -        fread(&d_temp,sizeof(double),1,fp);
  3837. -    else
  3838. -    {
  3839. -        fread(&f_temp,sizeof(float),1,fp);
  3840. -        d_temp = f_temp;
  3841. -    }
  3842. -    if ( o_flag == ROW_ORDER )
  3843. -        A->me[i / A->n][i % A->n].re = d_temp;
  3844. -    else if ( o_flag == COL_ORDER )
  3845. -        A->me[i % A->m][i / A->m].re = d_temp;
  3846. -    else
  3847. -        error(E_FORMAT,"zm_load");
  3848. -    }
  3849. -    
  3850. -    if ( mat.imag )         /* skip imaginary part */
  3851. -    for ( i = 0; i < A->m*A->n; i++ )
  3852. -    {
  3853. -        if ( p_flag == DOUBLE_PREC )
  3854. -        fread(&d_temp,sizeof(double),1,fp);
  3855. -        else
  3856. -        {
  3857. -        fread(&f_temp,sizeof(float),1,fp);
  3858. -        d_temp = f_temp;
  3859. -        }
  3860. -        if ( o_flag == ROW_ORDER )
  3861. -        A->me[i / A->n][i % A->n].im = d_temp;
  3862. -        else if ( o_flag == COL_ORDER )
  3863. -        A->me[i % A->m][i / A->m].im = d_temp;
  3864. -        else
  3865. -        error(E_FORMAT,"zm_load");
  3866. -    }
  3867. -    
  3868. -    return A;
  3869. -}
  3870. -
  3871. //GO.SYSIN DD zmatlab.c
  3872. echo zhsehldr.c 1>&2
  3873. sed >zhsehldr.c <<'//GO.SYSIN DD zhsehldr.c' 's/^-//'
  3874. -
  3875. -/**************************************************************************
  3876. -**
  3877. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  3878. -**
  3879. -**                 Meschach Library
  3880. -** 
  3881. -** This Meschach Library is provided "as is" without any express 
  3882. -** or implied warranty of any kind with respect to this software. 
  3883. -** In particular the authors shall not be liable for any direct, 
  3884. -** indirect, special, incidental or consequential damages arising 
  3885. -** in any way from use of the software.
  3886. -** 
  3887. -** Everyone is granted permission to copy, modify and redistribute this
  3888. -** Meschach Library, provided:
  3889. -**  1.  All copies contain this copyright notice.
  3890. -**  2.  All modified copies shall carry a notice stating who
  3891. -**      made the last modification and the date of such modification.
  3892. -**  3.  No charge is made for this software or works derived from it.  
  3893. -**      This clause shall not be construed as constraining other software
  3894. -**      distributed on the same medium as this software, nor is a
  3895. -**      distribution fee considered a charge.
  3896. -**
  3897. -***************************************************************************/
  3898. -
  3899. -
  3900. -/*
  3901. -        Files for matrix computations
  3902. -
  3903. -    Householder transformation file. Contains routines for calculating
  3904. -    householder transformations, applying them to vectors and matrices
  3905. -    by both row & column.
  3906. -
  3907. -    Complex version
  3908. -*/
  3909. -
  3910. -static    char    rcsid[] = "$Id: zhsehldr.c,v 1.2 1994/04/07 01:43:47 des Exp $";
  3911. -
  3912. -#include    <stdio.h>
  3913. -#include    <math.h>
  3914. -#include    "zmatrix.h"
  3915. -#include        "zmatrix2.h"
  3916. -
  3917. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  3918. -
  3919. -/* zhhvec -- calulates Householder vector to eliminate all entries after the
  3920. -    i0 entry of the vector vec. It is returned as out. May be in-situ */
  3921. -ZVEC    *zhhvec(vec,i0,beta,out,newval)
  3922. -ZVEC    *vec,*out;
  3923. -int    i0;
  3924. -Real    *beta;
  3925. -complex    *newval;
  3926. -{
  3927. -    complex    tmp;
  3928. -    Real    norm, abs_val;
  3929. -
  3930. -    if ( i0 < 0 || i0 >= vec->dim )
  3931. -        error(E_BOUNDS,"zhhvec");
  3932. -    out = _zv_copy(vec,out,i0);
  3933. -    tmp = _zin_prod(out,out,i0,Z_CONJ);
  3934. -    if ( tmp.re <= 0.0 )
  3935. -    {
  3936. -        *beta = 0.0;
  3937. -        *newval = out->ve[i0];
  3938. -        return (out);
  3939. -    }
  3940. -    norm = sqrt(tmp.re);
  3941. -    abs_val = zabs(out->ve[i0]);
  3942. -    *beta = 1.0/(norm * (norm+abs_val));
  3943. -    if ( abs_val == 0.0 )
  3944. -    {
  3945. -      newval->re = norm;
  3946. -      newval->im = 0.0;
  3947. -    }
  3948. -    else
  3949. -    { 
  3950. -      abs_val = -norm / abs_val;
  3951. -      newval->re = abs_val*out->ve[i0].re;
  3952. -      newval->im = abs_val*out->ve[i0].im;
  3953. -    }    abs_val = -norm / abs_val;
  3954. -    out->ve[i0].re -= newval->re;
  3955. -    out->ve[i0].im -= newval->im;
  3956. -
  3957. -    return (out);
  3958. -}
  3959. -
  3960. -/* zhhtrvec -- apply Householder transformation to vector -- may be in-situ */
  3961. -ZVEC    *zhhtrvec(hh,beta,i0,in,out)
  3962. -ZVEC    *hh,*in,*out;    /* hh = Householder vector */
  3963. -int    i0;
  3964. -double    beta;
  3965. -{
  3966. -    complex    scale, tmp;
  3967. -    /* u_int    i; */
  3968. -
  3969. -    if ( hh==ZVNULL || in==ZVNULL )
  3970. -        error(E_NULL,"zhhtrvec");
  3971. -    if ( in->dim != hh->dim )
  3972. -        error(E_SIZES,"zhhtrvec");
  3973. -    if ( i0 < 0 || i0 > in->dim )
  3974. -        error(E_BOUNDS,"zhhvec");
  3975. -
  3976. -    tmp = _zin_prod(hh,in,i0,Z_CONJ);
  3977. -    scale.re = -beta*tmp.re;
  3978. -    scale.im = -beta*tmp.im;
  3979. -    out = zv_copy(in,out);
  3980. -    __zmltadd__(&(out->ve[i0]),&(hh->ve[i0]),scale,
  3981. -            (int)(in->dim-i0),Z_NOCONJ);
  3982. -    /************************************************************
  3983. -    for ( i=i0; i<in->dim; i++ )
  3984. -        out->ve[i] = in->ve[i] - scale*hh->ve[i];
  3985. -    ************************************************************/
  3986. -
  3987. -    return (out);
  3988. -}
  3989. -
  3990. -/* zhhtrrows -- transform a matrix by a Householder vector by rows
  3991. -    starting at row i0 from column j0 -- in-situ */
  3992. -ZMAT    *zhhtrrows(M,i0,j0,hh,beta)
  3993. -ZMAT    *M;
  3994. -int    i0, j0;
  3995. -ZVEC    *hh;
  3996. -double    beta;
  3997. -{
  3998. -    complex    ip, scale;
  3999. -    int    i /*, j */;
  4000. -
  4001. -    if ( M==ZMNULL || hh==ZVNULL )
  4002. -        error(E_NULL,"zhhtrrows");
  4003. -    if ( M->n != hh->dim )
  4004. -        error(E_RANGE,"zhhtrrows");
  4005. -    if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  4006. -        error(E_BOUNDS,"zhhtrrows");
  4007. -
  4008. -    if ( beta == 0.0 )    return (M);
  4009. -
  4010. -    /* for each row ... */
  4011. -    for ( i = i0; i < M->m; i++ )
  4012. -    {    /* compute inner product */
  4013. -        ip = __zip__(&(M->me[i][j0]),&(hh->ve[j0]),
  4014. -                 (int)(M->n-j0),Z_NOCONJ);
  4015. -        /**************************************************
  4016. -        ip = 0.0;
  4017. -        for ( j = j0; j < M->n; j++ )
  4018. -            ip += M->me[i][j]*hh->ve[j];
  4019. -        **************************************************/
  4020. -        scale.re = -beta*ip.re;
  4021. -        scale.im = -beta*ip.im;
  4022. -        /* if ( scale == 0.0 ) */
  4023. -        if ( is_zero(scale) )
  4024. -            continue;
  4025. -
  4026. -        /* do operation */
  4027. -        __zmltadd__(&(M->me[i][j0]),&(hh->ve[j0]),scale,
  4028. -                (int)(M->n-j0),Z_CONJ);
  4029. -        /**************************************************
  4030. -        for ( j = j0; j < M->n; j++ )
  4031. -            M->me[i][j] -= scale*hh->ve[j];
  4032. -        **************************************************/
  4033. -    }
  4034. -
  4035. -    return (M);
  4036. -}
  4037. -
  4038. -
  4039. -/* zhhtrcols -- transform a matrix by a Householder vector by columns
  4040. -    starting at row i0 from column j0 -- in-situ */
  4041. -ZMAT    *zhhtrcols(M,i0,j0,hh,beta)
  4042. -ZMAT    *M;
  4043. -int    i0, j0;
  4044. -ZVEC    *hh;
  4045. -double    beta;
  4046. -{
  4047. -    /* Real    ip, scale; */
  4048. -    complex    scale;
  4049. -    int    i /*, k */;
  4050. -    static    ZVEC    *w = ZVNULL;
  4051. -
  4052. -    if ( M==ZMNULL || hh==ZVNULL )
  4053. -        error(E_NULL,"zhhtrcols");
  4054. -    if ( M->m != hh->dim )
  4055. -        error(E_SIZES,"zhhtrcols");
  4056. -    if ( i0 < 0 || i0 > M->m || j0 < 0 || j0 > M->n )
  4057. -        error(E_BOUNDS,"zhhtrcols");
  4058. -
  4059. -    if ( beta == 0.0 )    return (M);
  4060. -
  4061. -    w = zv_resize(w,M->n);
  4062. -    MEM_STAT_REG(w,TYPE_ZVEC);
  4063. -    zv_zero(w);
  4064. -
  4065. -    for ( i = i0; i < M->m; i++ )
  4066. -        /* if ( hh->ve[i] != 0.0 ) */
  4067. -        if ( ! is_zero(hh->ve[i]) )
  4068. -        __zmltadd__(&(w->ve[j0]),&(M->me[i][j0]),hh->ve[i],
  4069. -                (int)(M->n-j0),Z_CONJ);
  4070. -    for ( i = i0; i < M->m; i++ )
  4071. -        /* if ( hh->ve[i] != 0.0 ) */
  4072. -        if ( ! is_zero(hh->ve[i]) )
  4073. -        {
  4074. -        scale.re = -beta*hh->ve[i].re;
  4075. -        scale.im = -beta*hh->ve[i].im;
  4076. -        __zmltadd__(&(M->me[i][j0]),&(w->ve[j0]),scale,
  4077. -                (int)(M->n-j0),Z_CONJ);
  4078. -        }
  4079. -    return (M);
  4080. -}
  4081. -
  4082. //GO.SYSIN DD zhsehldr.c
  4083. echo zqrfctr.c 1>&2
  4084. sed >zqrfctr.c <<'//GO.SYSIN DD zqrfctr.c' 's/^-//'
  4085. -
  4086. -/**************************************************************************
  4087. -**
  4088. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4089. -**
  4090. -**                 Meschach Library
  4091. -** 
  4092. -** This Meschach Library is provided "as is" without any express 
  4093. -** or implied warranty of any kind with respect to this software. 
  4094. -** In particular the authors shall not be liable for any direct, 
  4095. -** indirect, special, incidental or consequential damages arising 
  4096. -** in any way from use of the software.
  4097. -** 
  4098. -** Everyone is granted permission to copy, modify and redistribute this
  4099. -** Meschach Library, provided:
  4100. -**  1.  All copies contain this copyright notice.
  4101. -**  2.  All modified copies shall carry a notice stating who
  4102. -**      made the last modification and the date of such modification.
  4103. -**  3.  No charge is made for this software or works derived from it.  
  4104. -**      This clause shall not be construed as constraining other software
  4105. -**      distributed on the same medium as this software, nor is a
  4106. -**      distribution fee considered a charge.
  4107. -**
  4108. -***************************************************************************/
  4109. -
  4110. -/*
  4111. -  This file contains the routines needed to perform QR factorisation
  4112. -  of matrices, as well as Householder transformations.
  4113. -  The internal "factored form" of a matrix A is not quite standard.
  4114. -  The diagonal of A is replaced by the diagonal of R -- not by the 1st non-zero
  4115. -  entries of the Householder vectors. The 1st non-zero entries are held in
  4116. -  the diag parameter of QRfactor(). The reason for this non-standard
  4117. -  representation is that it enables direct use of the Usolve() function
  4118. -  rather than requiring that  a seperate function be written just for this case.
  4119. -  See, e.g., QRsolve() below for more details.
  4120. -
  4121. -  Complex version
  4122. -  
  4123. -*/
  4124. -
  4125. -static    char    rcsid[] = "$Id: zqrfctr.c,v 1.1 1994/01/13 04:21:22 des Exp $";
  4126. -
  4127. -#include    <stdio.h>
  4128. -#include    <math.h>
  4129. -#include    "zmatrix.h"
  4130. -#include    "zmatrix2.h" 
  4131. -
  4132. -
  4133. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  4134. -
  4135. -
  4136. -#define        sign(x)    ((x) > 0.0 ? 1 : ((x) < 0.0 ? -1 : 0 ))
  4137. -
  4138. -/* Note: The usual representation of a Householder transformation is taken
  4139. -   to be:
  4140. -   P = I - beta.u.u*
  4141. -   where beta = 2/(u*.u) and u is called the Householder vector
  4142. -   (u* is the conjugate transposed vector of u
  4143. -*/
  4144. -
  4145. -/* zQRfactor -- forms the QR factorisation of A
  4146. -    -- factorisation stored in compact form as described above
  4147. -    (not quite standard format) */
  4148. -ZMAT    *zQRfactor(A,diag)
  4149. -ZMAT    *A;
  4150. -ZVEC    *diag;
  4151. -{
  4152. -    u_int    k,limit;
  4153. -    Real    beta;
  4154. -    static    ZVEC    *tmp1=ZVNULL;
  4155. -    
  4156. -    if ( ! A || ! diag )
  4157. -    error(E_NULL,"zQRfactor");
  4158. -    limit = min(A->m,A->n);
  4159. -    if ( diag->dim < limit )
  4160. -    error(E_SIZES,"zQRfactor");
  4161. -    
  4162. -    tmp1 = zv_resize(tmp1,A->m);
  4163. -    MEM_STAT_REG(tmp1,TYPE_ZVEC);
  4164. -    
  4165. -    for ( k=0; k<limit; k++ )
  4166. -    {
  4167. -    /* get H/holder vector for the k-th column */
  4168. -    zget_col(A,k,tmp1);
  4169. -    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
  4170. -    zhhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
  4171. -    diag->ve[k] = tmp1->ve[k];
  4172. -    
  4173. -    /* apply H/holder vector to remaining columns */
  4174. -    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
  4175. -    tracecatch(zhhtrcols(A,k,k+1,tmp1,beta),"zQRfactor");
  4176. -    }
  4177. -
  4178. -    return (A);
  4179. -}
  4180. -
  4181. -/* zQRCPfactor -- forms the QR factorisation of A with column pivoting
  4182. -   -- factorisation stored in compact form as described above
  4183. -   ( not quite standard format )                */
  4184. -ZMAT    *zQRCPfactor(A,diag,px)
  4185. -ZMAT    *A;
  4186. -ZVEC    *diag;
  4187. -PERM    *px;
  4188. -{
  4189. -    u_int    i, i_max, j, k, limit;
  4190. -    static    ZVEC    *tmp1=ZVNULL, *tmp2=ZVNULL;
  4191. -    static    VEC    *gamma=VNULL;
  4192. -    Real     beta;
  4193. -    Real    maxgamma, sum, tmp;
  4194. -    complex    ztmp;
  4195. -    
  4196. -    if ( ! A || ! diag || ! px )
  4197. -    error(E_NULL,"QRCPfactor");
  4198. -    limit = min(A->m,A->n);
  4199. -    if ( diag->dim < limit || px->size != A->n )
  4200. -    error(E_SIZES,"QRCPfactor");
  4201. -    
  4202. -    tmp1 = zv_resize(tmp1,A->m);
  4203. -    tmp2 = zv_resize(tmp2,A->m);
  4204. -    gamma = v_resize(gamma,A->n);
  4205. -    MEM_STAT_REG(tmp1,TYPE_ZVEC);
  4206. -    MEM_STAT_REG(tmp2,TYPE_ZVEC);
  4207. -    MEM_STAT_REG(gamma,TYPE_VEC);
  4208. -    
  4209. -    /* initialise gamma and px */
  4210. -    for ( j=0; j<A->n; j++ )
  4211. -    {
  4212. -    px->pe[j] = j;
  4213. -    sum = 0.0;
  4214. -    for ( i=0; i<A->m; i++ )
  4215. -        sum += square(A->me[i][j].re) + square(A->me[i][j].im);
  4216. -    gamma->ve[j] = sum;
  4217. -    }
  4218. -    
  4219. -    for ( k=0; k<limit; k++ )
  4220. -    {
  4221. -    /* find "best" column to use */
  4222. -    i_max = k;    maxgamma = gamma->ve[k];
  4223. -    for ( i=k+1; i<A->n; i++ )
  4224. -        /* Loop invariant:maxgamma=gamma[i_max]
  4225. -           >=gamma[l];l=k,...,i-1 */
  4226. -        if ( gamma->ve[i] > maxgamma )
  4227. -        {    maxgamma = gamma->ve[i]; i_max = i;    }
  4228. -    
  4229. -    /* swap columns if necessary */
  4230. -    if ( i_max != k )
  4231. -    {
  4232. -        /* swap gamma values */
  4233. -        tmp = gamma->ve[k];
  4234. -        gamma->ve[k] = gamma->ve[i_max];
  4235. -        gamma->ve[i_max] = tmp;
  4236. -        
  4237. -        /* update column permutation */
  4238. -        px_transp(px,k,i_max);
  4239. -        
  4240. -        /* swap columns of A */
  4241. -        for ( i=0; i<A->m; i++ )
  4242. -        {
  4243. -        ztmp = A->me[i][k];
  4244. -        A->me[i][k] = A->me[i][i_max];
  4245. -        A->me[i][i_max] = ztmp;
  4246. -        }
  4247. -    }
  4248. -    
  4249. -    /* get H/holder vector for the k-th column */
  4250. -    zget_col(A,k,tmp1);
  4251. -    /* hhvec(tmp1,k,&beta->ve[k],tmp1,&A->me[k][k]); */
  4252. -    zhhvec(tmp1,k,&beta,tmp1,&A->me[k][k]);
  4253. -    diag->ve[k] = tmp1->ve[k];
  4254. -    
  4255. -    /* apply H/holder vector to remaining columns */
  4256. -    /* hhtrcols(A,k,k+1,tmp1,beta->ve[k]); */
  4257. -    zhhtrcols(A,k,k+1,tmp1,beta);
  4258. -    
  4259. -    /* update gamma values */
  4260. -    for ( j=k+1; j<A->n; j++ )
  4261. -        gamma->ve[j] -= square(A->me[k][j].re)+square(A->me[k][j].im);
  4262. -    }
  4263. -
  4264. -    return (A);
  4265. -}
  4266. -
  4267. -/* zQsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact
  4268. -    form a la QRfactor()
  4269. -    -- may be in-situ */
  4270. -ZVEC    *_zQsolve(QR,diag,b,x,tmp)
  4271. -ZMAT    *QR;
  4272. -ZVEC    *diag, *b, *x, *tmp;
  4273. -{
  4274. -    u_int    dynamic;
  4275. -    int        k, limit;
  4276. -    Real    beta, r_ii, tmp_val;
  4277. -    
  4278. -    limit = min(QR->m,QR->n);
  4279. -    dynamic = FALSE;
  4280. -    if ( ! QR || ! diag || ! b )
  4281. -    error(E_NULL,"_zQsolve");
  4282. -    if ( diag->dim < limit || b->dim != QR->m )
  4283. -    error(E_SIZES,"_zQsolve");
  4284. -    x = zv_resize(x,QR->m);
  4285. -    if ( tmp == ZVNULL )
  4286. -    dynamic = TRUE;
  4287. -    tmp = zv_resize(tmp,QR->m);
  4288. -    
  4289. -    /* apply H/holder transforms in normal order */
  4290. -    x = zv_copy(b,x);
  4291. -    for ( k = 0 ; k < limit ; k++ )
  4292. -    {
  4293. -    zget_col(QR,k,tmp);
  4294. -    r_ii = zabs(tmp->ve[k]);
  4295. -    tmp->ve[k] = diag->ve[k];
  4296. -    tmp_val = (r_ii*zabs(diag->ve[k]));
  4297. -    beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  4298. -    /* hhtrvec(tmp,beta->ve[k],k,x,x); */
  4299. -    zhhtrvec(tmp,beta,k,x,x);
  4300. -    }
  4301. -    
  4302. -    if ( dynamic )
  4303. -    ZV_FREE(tmp);
  4304. -    
  4305. -    return (x);
  4306. -}
  4307. -
  4308. -/* zmakeQ -- constructs orthogonal matrix from Householder vectors stored in
  4309. -   compact QR form */
  4310. -ZMAT    *zmakeQ(QR,diag,Qout)
  4311. -ZMAT    *QR,*Qout;
  4312. -ZVEC    *diag;
  4313. -{
  4314. -    static    ZVEC    *tmp1=ZVNULL,*tmp2=ZVNULL;
  4315. -    u_int    i, limit;
  4316. -    Real    beta, r_ii, tmp_val;
  4317. -    int    j;
  4318. -
  4319. -    limit = min(QR->m,QR->n);
  4320. -    if ( ! QR || ! diag )
  4321. -    error(E_NULL,"zmakeQ");
  4322. -    if ( diag->dim < limit )
  4323. -    error(E_SIZES,"zmakeQ");
  4324. -    Qout = zm_resize(Qout,QR->m,QR->m);
  4325. -
  4326. -    tmp1 = zv_resize(tmp1,QR->m);    /* contains basis vec & columns of Q */
  4327. -    tmp2 = zv_resize(tmp2,QR->m);    /* contains H/holder vectors */
  4328. -    MEM_STAT_REG(tmp1,TYPE_ZVEC);
  4329. -    MEM_STAT_REG(tmp2,TYPE_ZVEC);
  4330. -
  4331. -    for ( i=0; i<QR->m ; i++ )
  4332. -    {    /* get i-th column of Q */
  4333. -    /* set up tmp1 as i-th basis vector */
  4334. -    for ( j=0; j<QR->m ; j++ )
  4335. -        tmp1->ve[j].re = tmp1->ve[j].im = 0.0;
  4336. -    tmp1->ve[i].re = 1.0;
  4337. -    
  4338. -    /* apply H/h transforms in reverse order */
  4339. -    for ( j=limit-1; j>=0; j-- )
  4340. -    {
  4341. -        zget_col(QR,j,tmp2);
  4342. -        r_ii = zabs(tmp2->ve[j]);
  4343. -        tmp2->ve[j] = diag->ve[j];
  4344. -        tmp_val = (r_ii*zabs(diag->ve[j]));
  4345. -        beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  4346. -        /* hhtrvec(tmp2,beta->ve[j],j,tmp1,tmp1); */
  4347. -        zhhtrvec(tmp2,beta,j,tmp1,tmp1);
  4348. -    }
  4349. -    
  4350. -    /* insert into Q */
  4351. -    zset_col(Qout,i,tmp1);
  4352. -    }
  4353. -
  4354. -    return (Qout);
  4355. -}
  4356. -
  4357. -/* zmakeR -- constructs upper triangular matrix from QR (compact form)
  4358. -    -- may be in-situ (all it does is zero the lower 1/2) */
  4359. -ZMAT    *zmakeR(QR,Rout)
  4360. -ZMAT    *QR,*Rout;
  4361. -{
  4362. -    u_int    i,j;
  4363. -    
  4364. -    if ( QR==ZMNULL )
  4365. -    error(E_NULL,"zmakeR");
  4366. -    Rout = zm_copy(QR,Rout);
  4367. -    
  4368. -    for ( i=1; i<QR->m; i++ )
  4369. -    for ( j=0; j<QR->n && j<i; j++ )
  4370. -        Rout->me[i][j].re = Rout->me[i][j].im = 0.0;
  4371. -    
  4372. -    return (Rout);
  4373. -}
  4374. -
  4375. -/* zQRsolve -- solves the system Q.R.x=b where Q & R are stored in compact form
  4376. -   -- returns x, which is created if necessary */
  4377. -ZVEC    *zQRsolve(QR,diag,b,x)
  4378. -ZMAT    *QR;
  4379. -ZVEC    *diag, *b, *x;
  4380. -{
  4381. -    int    limit;
  4382. -    static    ZVEC    *tmp = ZVNULL;
  4383. -    
  4384. -    if ( ! QR || ! diag || ! b )
  4385. -    error(E_NULL,"zQRsolve");
  4386. -    limit = min(QR->m,QR->n);
  4387. -    if ( diag->dim < limit || b->dim != QR->m )
  4388. -    error(E_SIZES,"zQRsolve");
  4389. -    tmp = zv_resize(tmp,limit);
  4390. -    MEM_STAT_REG(tmp,TYPE_ZVEC);
  4391. -
  4392. -    x = zv_resize(x,QR->n);
  4393. -    _zQsolve(QR,diag,b,x,tmp);
  4394. -    x = zUsolve(QR,x,x,0.0);
  4395. -    x = zv_resize(x,QR->n);
  4396. -
  4397. -    return x;
  4398. -}
  4399. -
  4400. -/* zQRAsolve -- solves the system (Q.R)*.x = b
  4401. -    -- Q & R are stored in compact form
  4402. -    -- returns x, which is created if necessary */
  4403. -ZVEC    *zQRAsolve(QR,diag,b,x)
  4404. -ZMAT    *QR;
  4405. -ZVEC    *diag, *b, *x;
  4406. -{
  4407. -    int        j, limit;
  4408. -    Real    beta, r_ii, tmp_val;
  4409. -    static    ZVEC    *tmp = ZVNULL;
  4410. -    
  4411. -    if ( ! QR || ! diag || ! b )
  4412. -    error(E_NULL,"zQRAsolve");
  4413. -    limit = min(QR->m,QR->n);
  4414. -    if ( diag->dim < limit || b->dim != QR->n )
  4415. -    error(E_SIZES,"zQRAsolve");
  4416. -
  4417. -    x = zv_resize(x,QR->m);
  4418. -    x = zUAsolve(QR,b,x,0.0);
  4419. -    x = zv_resize(x,QR->m);
  4420. -
  4421. -    tmp = zv_resize(tmp,x->dim);
  4422. -    MEM_STAT_REG(tmp,TYPE_ZVEC);
  4423. -    printf("zQRAsolve: tmp->dim = %d, x->dim = %d\n", tmp->dim, x->dim);
  4424. -    
  4425. -    /* apply H/h transforms in reverse order */
  4426. -    for ( j=limit-1; j>=0; j-- )
  4427. -    {
  4428. -    zget_col(QR,j,tmp);
  4429. -    tmp = zv_resize(tmp,QR->m);
  4430. -    r_ii = zabs(tmp->ve[j]);
  4431. -    tmp->ve[j] = diag->ve[j];
  4432. -    tmp_val = (r_ii*zabs(diag->ve[j]));
  4433. -    beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  4434. -    zhhtrvec(tmp,beta,j,x,x);
  4435. -    }
  4436. -
  4437. -
  4438. -    return x;
  4439. -}
  4440. -
  4441. -/* zQRCPsolve -- solves A.x = b where A is factored by QRCPfactor()
  4442. -   -- assumes that A is in the compact factored form */
  4443. -ZVEC    *zQRCPsolve(QR,diag,pivot,b,x)
  4444. -ZMAT    *QR;
  4445. -ZVEC    *diag;
  4446. -PERM    *pivot;
  4447. -ZVEC    *b, *x;
  4448. -{
  4449. -    if ( ! QR || ! diag || ! pivot || ! b )
  4450. -    error(E_NULL,"zQRCPsolve");
  4451. -    if ( (QR->m > diag->dim && QR->n > diag->dim) || QR->n != pivot->size )
  4452. -    error(E_SIZES,"zQRCPsolve");
  4453. -    
  4454. -    x = zQRsolve(QR,diag,b,x);
  4455. -    x = pxinv_zvec(pivot,x,x);
  4456. -
  4457. -    return x;
  4458. -}
  4459. -
  4460. -/* zUmlt -- compute out = upper_triang(U).x
  4461. -    -- may be in situ */
  4462. -ZVEC    *zUmlt(U,x,out)
  4463. -ZMAT    *U;
  4464. -ZVEC    *x, *out;
  4465. -{
  4466. -    int        i, limit;
  4467. -
  4468. -    if ( U == ZMNULL || x == ZVNULL )
  4469. -    error(E_NULL,"zUmlt");
  4470. -    limit = min(U->m,U->n);
  4471. -    if ( limit != x->dim )
  4472. -    error(E_SIZES,"zUmlt");
  4473. -    if ( out == ZVNULL || out->dim < limit )
  4474. -    out = zv_resize(out,limit);
  4475. -
  4476. -    for ( i = 0; i < limit; i++ )
  4477. -    out->ve[i] = __zip__(&(x->ve[i]),&(U->me[i][i]),limit - i,Z_NOCONJ);
  4478. -    return out;
  4479. -}
  4480. -
  4481. -/* zUAmlt -- returns out = upper_triang(U)^T.x */
  4482. -ZVEC    *zUAmlt(U,x,out)
  4483. -ZMAT    *U;
  4484. -ZVEC    *x, *out;
  4485. -{
  4486. -    /* complex    sum; */
  4487. -    complex    tmp;
  4488. -    int        i, limit;
  4489. -
  4490. -    if ( U == ZMNULL || x == ZVNULL )
  4491. -    error(E_NULL,"zUAmlt");
  4492. -    limit = min(U->m,U->n);
  4493. -    if ( out == ZVNULL || out->dim < limit )
  4494. -    out = zv_resize(out,limit);
  4495. -
  4496. -    for ( i = limit-1; i >= 0; i-- )
  4497. -    {
  4498. -    tmp = x->ve[i];
  4499. -    out->ve[i].re = out->ve[i].im = 0.0;
  4500. -    __zmltadd__(&(out->ve[i]),&(U->me[i][i]),tmp,limit-i-1,Z_CONJ);
  4501. -    }
  4502. -
  4503. -    return out;
  4504. -}
  4505. -
  4506. -
  4507. -/* zQRcondest -- returns an estimate of the 2-norm condition number of the
  4508. -        matrix factorised by QRfactor() or QRCPfactor()
  4509. -    -- note that as Q does not affect the 2-norm condition number,
  4510. -        it is not necessary to pass the diag, beta (or pivot) vectors
  4511. -    -- generates a lower bound on the true condition number
  4512. -    -- if the matrix is exactly singular, HUGE is returned
  4513. -    -- note that QRcondest() is likely to be more reliable for
  4514. -        matrices factored using QRCPfactor() */
  4515. -double    zQRcondest(QR)
  4516. -ZMAT    *QR;
  4517. -{
  4518. -    static    ZVEC    *y=ZVNULL;
  4519. -    Real    norm, norm1, norm2, tmp1, tmp2;
  4520. -    complex    sum, tmp;
  4521. -    int        i, j, limit;
  4522. -
  4523. -    if ( QR == ZMNULL )
  4524. -    error(E_NULL,"zQRcondest");
  4525. -
  4526. -    limit = min(QR->m,QR->n);
  4527. -    for ( i = 0; i < limit; i++ )
  4528. -    /* if ( QR->me[i][i] == 0.0 ) */
  4529. -    if ( is_zero(QR->me[i][i]) )
  4530. -        return HUGE;
  4531. -
  4532. -    y = zv_resize(y,limit);
  4533. -    MEM_STAT_REG(y,TYPE_ZVEC);
  4534. -    /* use the trick for getting a unit vector y with ||R.y||_inf small
  4535. -       from the LU condition estimator */
  4536. -    for ( i = 0; i < limit; i++ )
  4537. -    {
  4538. -    sum.re = sum.im = 0.0;
  4539. -    for ( j = 0; j < i; j++ )
  4540. -        /* sum -= QR->me[j][i]*y->ve[j]; */
  4541. -        sum = zsub(sum,zmlt(QR->me[j][i],y->ve[j]));
  4542. -    /* sum -= (sum < 0.0) ? 1.0 : -1.0; */
  4543. -    norm1 = zabs(sum);
  4544. -    if ( norm1 == 0.0 )
  4545. -        sum.re = 1.0;
  4546. -    else
  4547. -    {
  4548. -        sum.re += sum.re / norm1;
  4549. -        sum.im += sum.im / norm1;
  4550. -    }
  4551. -    /* y->ve[i] = sum / QR->me[i][i]; */
  4552. -    y->ve[i] = zdiv(sum,QR->me[i][i]);
  4553. -    }
  4554. -    zUAmlt(QR,y,y);
  4555. -
  4556. -    /* now apply inverse power method to R*.R */
  4557. -    for ( i = 0; i < 3; i++ )
  4558. -    {
  4559. -    tmp1 = zv_norm2(y);
  4560. -    zv_mlt(zmake(1.0/tmp1,0.0),y,y);
  4561. -    zUAsolve(QR,y,y,0.0);
  4562. -    tmp2 = zv_norm2(y);
  4563. -    zv_mlt(zmake(1.0/tmp2,0.0),y,y);
  4564. -    zUsolve(QR,y,y,0.0);
  4565. -    }
  4566. -    /* now compute approximation for ||R^{-1}||_2 */
  4567. -    norm1 = sqrt(tmp1)*sqrt(tmp2);
  4568. -
  4569. -    /* now use complementary approach to compute approximation to ||R||_2 */
  4570. -    for ( i = limit-1; i >= 0; i-- )
  4571. -    {
  4572. -    sum.re = sum.im = 0.0;
  4573. -    for ( j = i+1; j < limit; j++ )
  4574. -        sum = zadd(sum,zmlt(QR->me[i][j],y->ve[j]));
  4575. -    if ( is_zero(QR->me[i][i]) )
  4576. -        return HUGE;
  4577. -    tmp = zdiv(sum,QR->me[i][i]);
  4578. -    if ( is_zero(tmp) )
  4579. -    {
  4580. -        y->ve[i].re = 1.0;
  4581. -        y->ve[i].im = 0.0;
  4582. -    }
  4583. -    else
  4584. -    {
  4585. -        norm = zabs(tmp);
  4586. -        y->ve[i].re = sum.re / norm;
  4587. -        y->ve[i].im = sum.im / norm;
  4588. -    }
  4589. -    /* y->ve[i] = (sum >= 0.0) ? 1.0 : -1.0; */
  4590. -    /* y->ve[i] = (QR->me[i][i] >= 0.0) ? y->ve[i] : - y->ve[i]; */
  4591. -    }
  4592. -
  4593. -    /* now apply power method to R*.R */
  4594. -    for ( i = 0; i < 3; i++ )
  4595. -    {
  4596. -    tmp1 = zv_norm2(y);
  4597. -    zv_mlt(zmake(1.0/tmp1,0.0),y,y);
  4598. -    zUmlt(QR,y,y);
  4599. -    tmp2 = zv_norm2(y);
  4600. -    zv_mlt(zmake(1.0/tmp2,0.0),y,y);
  4601. -    zUAmlt(QR,y,y);
  4602. -    }
  4603. -    norm2 = sqrt(tmp1)*sqrt(tmp2);
  4604. -
  4605. -    /* printf("QRcondest: norm1 = %g, norm2 = %g\n",norm1,norm2); */
  4606. -
  4607. -    return norm1*norm2;
  4608. -}
  4609. -
  4610. //GO.SYSIN DD zqrfctr.c
  4611. echo zgivens.c 1>&2
  4612. sed >zgivens.c <<'//GO.SYSIN DD zgivens.c' 's/^-//'
  4613. -
  4614. -/**************************************************************************
  4615. -**
  4616. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4617. -**
  4618. -**                 Meschach Library
  4619. -** 
  4620. -** This Meschach Library is provided "as is" without any express 
  4621. -** or implied warranty of any kind with respect to this software. 
  4622. -** In particular the authors shall not be liable for any direct, 
  4623. -** indirect, special, incidental or consequential damages arising 
  4624. -** in any way from use of the software.
  4625. -** 
  4626. -** Everyone is granted permission to copy, modify and redistribute this
  4627. -** Meschach Library, provided:
  4628. -**  1.  All copies contain this copyright notice.
  4629. -**  2.  All modified copies shall carry a notice stating who
  4630. -**      made the last modification and the date of such modification.
  4631. -**  3.  No charge is made for this software or works derived from it.  
  4632. -**      This clause shall not be construed as constraining other software
  4633. -**      distributed on the same medium as this software, nor is a
  4634. -**      distribution fee considered a charge.
  4635. -**
  4636. -***************************************************************************/
  4637. -
  4638. -
  4639. -/*
  4640. -    Givens operations file. Contains routines for calculating and
  4641. -    applying givens rotations for/to vectors and also to matrices by
  4642. -    row and by column.
  4643. -
  4644. -    Complex version.
  4645. -*/
  4646. -
  4647. -static    char    rcsid[] = "$Id: ";
  4648. -
  4649. -#include    <stdio.h>
  4650. -#include    <math.h>
  4651. -#include    "zmatrix.h"
  4652. -#include        "zmatrix2.h"
  4653. -
  4654. -/*
  4655. -    (Complex) Givens rotation matrix:
  4656. -        [ c   -s ]
  4657. -        [ s*   c ]
  4658. -    Note that c is real and s is complex
  4659. -*/
  4660. -
  4661. -/* zgivens -- returns c,s parameters for Givens rotation to
  4662. -        eliminate y in the **column** vector [ x y ] */
  4663. -void    zgivens(x,y,c,s)
  4664. -complex    x,y,*s;
  4665. -Real    *c;
  4666. -{
  4667. -    Real    inv_norm, norm;
  4668. -    complex    tmp;
  4669. -
  4670. -    /* this is a safe way of computing sqrt(|x|^2+|y|^2) */
  4671. -    tmp.re = zabs(x);    tmp.im = zabs(y);
  4672. -    norm = zabs(tmp);
  4673. -
  4674. -    if ( norm == 0.0 )
  4675. -    {    *c = 1.0;    s->re = s->im = 0.0;    } /* identity */
  4676. -    else
  4677. -    {
  4678. -        inv_norm = 1.0 / tmp.re;    /* inv_norm = 1/|x| */
  4679. -        x.re *= inv_norm;
  4680. -        x.im *= inv_norm;        /* normalise x */
  4681. -        inv_norm = 1.0/norm;        /* inv_norm = 1/||[x,y]||2 */
  4682. -        *c = tmp.re * inv_norm;
  4683. -        /* now compute - conj(normalised x).y/||[x,y]||2 */
  4684. -        s->re = - inv_norm*(x.re*y.re + x.im*y.im);
  4685. -        s->im =   inv_norm*(x.re*y.im - x.im*y.re);
  4686. -    }
  4687. -}
  4688. -
  4689. -/* rot_zvec -- apply Givens rotation to x's i & k components */
  4690. -ZVEC    *rot_zvec(x,i,k,c,s,out)
  4691. -ZVEC    *x,*out;
  4692. -int    i,k;
  4693. -double    c;
  4694. -complex    s;
  4695. -{
  4696. -
  4697. -    complex    temp1, temp2;
  4698. -
  4699. -    if ( x==ZVNULL )
  4700. -        error(E_NULL,"rot_zvec");
  4701. -    if ( i < 0 || i >= x->dim || k < 0 || k >= x->dim )
  4702. -        error(E_RANGE,"rot_zvec");
  4703. -    if ( x != out )
  4704. -        out = zv_copy(x,out);
  4705. -
  4706. -    /* temp1 = c*out->ve[i] - s*out->ve[k]; */
  4707. -    temp1.re = c*out->ve[i].re
  4708. -        - s.re*out->ve[k].re + s.im*out->ve[k].im;
  4709. -    temp1.im = c*out->ve[i].im
  4710. -        - s.re*out->ve[k].im - s.im*out->ve[k].re;
  4711. -
  4712. -    /* temp2 = c*out->ve[k] + zconj(s)*out->ve[i]; */
  4713. -    temp2.re = c*out->ve[k].re
  4714. -        + s.re*out->ve[i].re + s.im*out->ve[i].im;
  4715. -    temp2.im = c*out->ve[k].im
  4716. -        + s.re*out->ve[i].im - s.im*out->ve[i].re;
  4717. -
  4718. -    out->ve[i] = temp1;
  4719. -    out->ve[k] = temp2;
  4720. -
  4721. -    return (out);
  4722. -}
  4723. -
  4724. -/* zrot_rows -- premultiply mat by givens rotation described by c,s */
  4725. -ZMAT    *zrot_rows(mat,i,k,c,s,out)
  4726. -ZMAT    *mat,*out;
  4727. -int    i,k;
  4728. -double    c;
  4729. -complex    s;
  4730. -{
  4731. -    u_int    j;
  4732. -    complex    temp1, temp2;
  4733. -
  4734. -    if ( mat==ZMNULL )
  4735. -        error(E_NULL,"zrot_rows");
  4736. -    if ( i < 0 || i >= mat->m || k < 0 || k >= mat->m )
  4737. -        error(E_RANGE,"zrot_rows");
  4738. -    out = zm_copy(mat,out);
  4739. -
  4740. -    /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  4741. -    for ( j=0; j<mat->n; j++ )
  4742. -    {
  4743. -        /* temp1 = c*out->me[i][j] - s*out->me[k][j]; */
  4744. -        temp1.re = c*out->me[i][j].re
  4745. -        - s.re*out->me[k][j].re + s.im*out->me[k][j].im;
  4746. -        temp1.im = c*out->me[i][j].im
  4747. -        - s.re*out->me[k][j].im - s.im*out->me[k][j].re;
  4748. -        
  4749. -        /* temp2 = c*out->me[k][j] + conj(s)*out->me[i][j]; */
  4750. -        temp2.re = c*out->me[k][j].re
  4751. -        + s.re*out->me[i][j].re + s.im*out->me[i][j].im;
  4752. -        temp2.im = c*out->me[k][j].im
  4753. -        + s.re*out->me[i][j].im - s.im*out->me[i][j].re;
  4754. -        
  4755. -        out->me[i][j] = temp1;
  4756. -        out->me[k][j] = temp2;
  4757. -    }
  4758. -
  4759. -    return (out);
  4760. -}
  4761. -
  4762. -/* zrot_cols -- postmultiply mat by adjoint Givens rotation described by c,s */
  4763. -ZMAT    *zrot_cols(mat,i,k,c,s,out)
  4764. -ZMAT    *mat,*out;
  4765. -int    i,k;
  4766. -double    c;
  4767. -complex    s;
  4768. -{
  4769. -    u_int    j;
  4770. -    complex    x, y;
  4771. -
  4772. -    if ( mat==ZMNULL )
  4773. -        error(E_NULL,"zrot_cols");
  4774. -    if ( i < 0 || i >= mat->n || k < 0 || k >= mat->n )
  4775. -        error(E_RANGE,"zrot_cols");
  4776. -    out = zm_copy(mat,out);
  4777. -
  4778. -    for ( j=0; j<mat->m; j++ )
  4779. -    {
  4780. -        x = out->me[j][i];    y = out->me[j][k];
  4781. -        /* out->me[j][i] = c*x - conj(s)*y; */
  4782. -        out->me[j][i].re = c*x.re - s.re*y.re - s.im*y.im;
  4783. -        out->me[j][i].im = c*x.im - s.re*y.im + s.im*y.re;
  4784. -        
  4785. -        /* out->me[j][k] = c*y + s*x; */
  4786. -        out->me[j][k].re = c*y.re + s.re*x.re - s.im*x.im;
  4787. -        out->me[j][k].im = c*y.im + s.re*x.im + s.im*x.re;
  4788. -    }
  4789. -
  4790. -    return (out);
  4791. -}
  4792. -
  4793. //GO.SYSIN DD zgivens.c
  4794. echo zhessen.c 1>&2
  4795. sed >zhessen.c <<'//GO.SYSIN DD zhessen.c' 's/^-//'
  4796. -
  4797. -/**************************************************************************
  4798. -**
  4799. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4800. -**
  4801. -**                 Meschach Library
  4802. -** 
  4803. -** This Meschach Library is provided "as is" without any express 
  4804. -** or implied warranty of any kind with respect to this software. 
  4805. -** In particular the authors shall not be liable for any direct, 
  4806. -** indirect, special, incidental or consequential damages arising 
  4807. -** in any way from use of the software.
  4808. -** 
  4809. -** Everyone is granted permission to copy, modify and redistribute this
  4810. -** Meschach Library, provided:
  4811. -**  1.  All copies contain this copyright notice.
  4812. -**  2.  All modified copies shall carry a notice stating who
  4813. -**      made the last modification and the date of such modification.
  4814. -**  3.  No charge is made for this software or works derived from it.  
  4815. -**      This clause shall not be construed as constraining other software
  4816. -**      distributed on the same medium as this software, nor is a
  4817. -**      distribution fee considered a charge.
  4818. -**
  4819. -***************************************************************************/
  4820. -
  4821. -
  4822. -/*
  4823. -        File containing routines for determining Hessenberg
  4824. -    factorisations.
  4825. -
  4826. -    Complex version
  4827. -*/
  4828. -
  4829. -static    char    rcsid[] = "$Id: zhessen.c,v 1.1 1994/01/13 04:27:00 des Exp $";
  4830. -
  4831. -#include    <stdio.h>
  4832. -#include    "zmatrix.h"
  4833. -#include        "zmatrix2.h"
  4834. -
  4835. -
  4836. -/* zHfactor -- compute Hessenberg factorisation in compact form.
  4837. -    -- factorisation performed in situ
  4838. -    -- for details of the compact form see zQRfactor.c and zmatrix2.doc */
  4839. -ZMAT    *zHfactor(A, diag)
  4840. -ZMAT    *A;
  4841. -ZVEC    *diag;
  4842. -{
  4843. -    static    ZVEC    *tmp1 = ZVNULL;
  4844. -    Real    beta;
  4845. -    int    k, limit;
  4846. -
  4847. -    if ( ! A || ! diag )
  4848. -        error(E_NULL,"zHfactor");
  4849. -    if ( diag->dim < A->m - 1 )
  4850. -        error(E_SIZES,"zHfactor");
  4851. -    if ( A->m != A->n )
  4852. -        error(E_SQUARE,"zHfactor");
  4853. -    limit = A->m - 1;
  4854. -
  4855. -    tmp1 = zv_resize(tmp1,A->m);
  4856. -    MEM_STAT_REG(tmp1,TYPE_ZVEC);
  4857. -
  4858. -    for ( k = 0; k < limit; k++ )
  4859. -    {
  4860. -        zget_col(A,k,tmp1);
  4861. -        zhhvec(tmp1,k+1,&beta,tmp1,&A->me[k+1][k]);
  4862. -        diag->ve[k] = tmp1->ve[k+1];
  4863. -        /* printf("zHfactor: k = %d, beta = %g, tmp1 =\n",k,beta);
  4864. -        zv_output(tmp1); */
  4865. -        
  4866. -        zhhtrcols(A,k+1,k+1,tmp1,beta);
  4867. -        zhhtrrows(A,0  ,k+1,tmp1,beta);
  4868. -        /* printf("# at stage k = %d, A =\n",k);    zm_output(A); */
  4869. -    }
  4870. -
  4871. -    return (A);
  4872. -}
  4873. -
  4874. -/* zHQunpack -- unpack the compact representation of H and Q of a
  4875. -    Hessenberg factorisation
  4876. -    -- if either H or Q is NULL, then it is not unpacked
  4877. -    -- it can be in situ with HQ == H
  4878. -    -- returns HQ
  4879. -*/
  4880. -ZMAT    *zHQunpack(HQ,diag,Q,H)
  4881. -ZMAT    *HQ, *Q, *H;
  4882. -ZVEC    *diag;
  4883. -{
  4884. -    int    i, j, limit;
  4885. -    Real    beta, r_ii, tmp_val;
  4886. -    static    ZVEC    *tmp1 = ZVNULL, *tmp2 = ZVNULL;
  4887. -
  4888. -    if ( HQ==ZMNULL || diag==ZVNULL )
  4889. -        error(E_NULL,"zHQunpack");
  4890. -    if ( HQ == Q || H == Q )
  4891. -        error(E_INSITU,"zHQunpack");
  4892. -    limit = HQ->m - 1;
  4893. -    if ( diag->dim < limit )
  4894. -        error(E_SIZES,"zHQunpack");
  4895. -    if ( HQ->m != HQ->n )
  4896. -        error(E_SQUARE,"zHQunpack");
  4897. -
  4898. -
  4899. -    if ( Q != ZMNULL )
  4900. -    {
  4901. -        Q = zm_resize(Q,HQ->m,HQ->m);
  4902. -        tmp1 = zv_resize(tmp1,H->m);
  4903. -        tmp2 = zv_resize(tmp2,H->m);
  4904. -        MEM_STAT_REG(tmp1,TYPE_ZVEC);
  4905. -        MEM_STAT_REG(tmp2,TYPE_ZVEC);
  4906. -        
  4907. -        for ( i = 0; i < H->m; i++ )
  4908. -        {
  4909. -        /* tmp1 = i'th basis vector */
  4910. -        for ( j = 0; j < H->m; j++ )
  4911. -            tmp1->ve[j].re = tmp1->ve[j].im = 0.0;
  4912. -        tmp1->ve[i].re = 1.0;
  4913. -        
  4914. -        /* apply H/h transforms in reverse order */
  4915. -        for ( j = limit-1; j >= 0; j-- )
  4916. -        {
  4917. -            zget_col(HQ,j,tmp2);
  4918. -            r_ii = zabs(tmp2->ve[j+1]);
  4919. -            tmp2->ve[j+1] = diag->ve[j];
  4920. -            tmp_val = (r_ii*zabs(diag->ve[j]));
  4921. -            beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
  4922. -            /* printf("zHQunpack: j = %d, beta = %g, tmp2 =\n",
  4923. -               j,beta);
  4924. -            zv_output(tmp2); */
  4925. -            zhhtrvec(tmp2,beta,j+1,tmp1,tmp1);
  4926. -        }
  4927. -        
  4928. -        /* insert into Q */
  4929. -        zset_col(Q,i,tmp1);
  4930. -        }
  4931. -    }
  4932. -
  4933. -    if ( H != ZMNULL )
  4934. -    {
  4935. -        H = zm_copy(HQ,H);
  4936. -        
  4937. -        limit = H->m;
  4938. -        for ( i = 1; i < limit; i++ )
  4939. -        for ( j = 0; j < i-1; j++ )
  4940. -            H->me[i][j].re = H->me[i][j].im = 0.0;
  4941. -    }
  4942. -
  4943. -    return HQ;
  4944. -}
  4945. -
  4946. -
  4947. -
  4948. //GO.SYSIN DD zhessen.c
  4949. echo zschur.c 1>&2
  4950. sed >zschur.c <<'//GO.SYSIN DD zschur.c' 's/^-//'
  4951. -
  4952. -/**************************************************************************
  4953. -**
  4954. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4955. -**
  4956. -**                 Meschach Library
  4957. -** 
  4958. -** This Meschach Library is provided "as is" without any express 
  4959. -** or implied warranty of any kind with respect to this software. 
  4960. -** In particular the authors shall not be liable for any direct, 
  4961. -** indirect, special, incidental or consequential damages arising 
  4962. -** in any way from use of the software.
  4963. -** 
  4964. -** Everyone is granted permission to copy, modify and redistribute this
  4965. -** Meschach Library, provided:
  4966. -**  1.  All copies contain this copyright notice.
  4967. -**  2.  All modified copies shall carry a notice stating who
  4968. -**      made the last modification and the date of such modification.
  4969. -**  3.  No charge is made for this software or works derived from it.  
  4970. -**      This clause shall not be construed as constraining other software
  4971. -**      distributed on the same medium as this software, nor is a
  4972. -**      distribution fee considered a charge.
  4973. -**
  4974. -***************************************************************************/
  4975. -
  4976. -/*    
  4977. -    File containing routines for computing the Schur decomposition
  4978. -    of a complex non-symmetric matrix
  4979. -    See also: hessen.c
  4980. -    Complex version
  4981. -*/
  4982. -
  4983. -
  4984. -#include    <stdio.h>
  4985. -#include    <math.h>
  4986. -#include    "zmatrix.h"
  4987. -#include        "zmatrix2.h"
  4988. -
  4989. -
  4990. -#define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  4991. -#define    b2s(t_or_f)    ((t_or_f) ? "TRUE" : "FALSE")
  4992. -
  4993. -
  4994. -/* zschur -- computes the Schur decomposition of the matrix A in situ
  4995. -    -- optionally, gives Q matrix such that Q^T.A.Q is upper triangular
  4996. -    -- returns upper triangular Schur matrix */
  4997. -ZMAT    *zschur(A,Q)
  4998. -ZMAT    *A, *Q;
  4999. -{
  5000. -    int        i, j, iter, k, k_min, k_max, k_tmp, n, split;
  5001. -    Real    c;
  5002. -    complex    det, discrim, lambda, lambda0, lambda1, s, sum, ztmp;
  5003. -    complex    x, y;    /* for chasing algorithm */
  5004. -    complex    **A_me;
  5005. -    static    ZVEC    *diag=ZVNULL;
  5006. -    
  5007. -    if ( ! A )
  5008. -    error(E_NULL,"zschur");
  5009. -    if ( A->m != A->n || ( Q && Q->m != Q->n ) )
  5010. -    error(E_SQUARE,"zschur");
  5011. -    if ( Q != ZMNULL && Q->m != A->m )
  5012. -    error(E_SIZES,"zschur");
  5013. -    n = A->n;
  5014. -    diag = zv_resize(diag,A->n);
  5015. -    MEM_STAT_REG(diag,TYPE_ZVEC);
  5016. -    /* compute Hessenberg form */
  5017. -    zHfactor(A,diag);
  5018. -    
  5019. -    /* save Q if necessary, and make A explicitly Hessenberg */
  5020. -    zHQunpack(A,diag,Q,A);
  5021. -
  5022. -    k_min = 0;    A_me = A->me;
  5023. -
  5024. -    while ( k_min < n )
  5025. -    {
  5026. -    /* find k_max to suit:
  5027. -       submatrix k_min..k_max should be irreducible */
  5028. -    k_max = n-1;
  5029. -    for ( k = k_min; k < k_max; k++ )
  5030. -        if ( is_zero(A_me[k+1][k]) )
  5031. -        {    k_max = k;    break;    }
  5032. -
  5033. -    if ( k_max <= k_min )
  5034. -    {
  5035. -        k_min = k_max + 1;
  5036. -        continue;        /* outer loop */
  5037. -    }
  5038. -
  5039. -    /* now have r x r block with r >= 2:
  5040. -       apply Francis QR step until block splits */
  5041. -    split = FALSE;        iter = 0;
  5042. -    while ( ! split )
  5043. -    {
  5044. -        complex    a00, a01, a10, a11;
  5045. -        iter++;
  5046. -        
  5047. -        /* set up Wilkinson/Francis complex shift */
  5048. -        /* use the smallest eigenvalue of the bottom 2 x 2 submatrix */
  5049. -        k_tmp = k_max - 1;
  5050. -
  5051. -        a00 = A_me[k_tmp][k_tmp];
  5052. -        a01 = A_me[k_tmp][k_max];
  5053. -        a10 = A_me[k_max][k_tmp];
  5054. -        a11 = A_me[k_max][k_max];
  5055. -        ztmp.re = 0.5*(a00.re - a11.re);
  5056. -        ztmp.im = 0.5*(a00.im - a11.im);
  5057. -        discrim = zsqrt(zadd(zmlt(ztmp,ztmp),zmlt(a01,a10)));
  5058. -        sum.re  = 0.5*(a00.re + a11.re);
  5059. -        sum.im  = 0.5*(a00.im + a11.im);
  5060. -        lambda0 = zadd(sum,discrim);
  5061. -        lambda1 = zsub(sum,discrim);
  5062. -        det = zsub(zmlt(a00,a11),zmlt(a01,a10));
  5063. -        if ( zabs(lambda0) > zabs(lambda1) )
  5064. -        lambda = zdiv(det,lambda0);
  5065. -        else
  5066. -        lambda = zdiv(det,lambda1);
  5067. -
  5068. -        /* perturb shift if convergence is slow */
  5069. -        if ( (iter % 10) == 0 )
  5070. -        {
  5071. -        lambda.re += iter*0.02;
  5072. -        lambda.im += iter*0.02;
  5073. -        }
  5074. -
  5075. -        /* set up Householder transformations */
  5076. -        k_tmp = k_min + 1;
  5077. -
  5078. -        x = zsub(A->me[k_min][k_min],lambda);
  5079. -        y = A->me[k_min+1][k_min];
  5080. -
  5081. -        /* use Givens' rotations to "chase" off-Hessenberg entry */
  5082. -        for ( k = k_min; k <= k_max-1; k++ )
  5083. -        {
  5084. -        zgivens(x,y,&c,&s);
  5085. -        zrot_cols(A,k,k+1,c,s,A);
  5086. -        zrot_rows(A,k,k+1,c,s,A);
  5087. -        if ( Q != ZMNULL )
  5088. -            zrot_cols(Q,k,k+1,c,s,Q);
  5089. -
  5090. -        /* zero things that should be zero */
  5091. -        if ( k > k_min )
  5092. -            A->me[k+1][k-1].re = A->me[k+1][k-1].im = 0.0;
  5093. -
  5094. -        /* get next entry to chase along sub-diagonal */
  5095. -        x = A->me[k+1][k];
  5096. -        if ( k <= k_max - 2 )
  5097. -            y = A->me[k+2][k];
  5098. -        else
  5099. -            y.re = y.im = 0.0;
  5100. -        }
  5101. -
  5102. -        for ( k = k_min; k <= k_max-2; k++ )
  5103. -        {
  5104. -        /* zero appropriate sub-diagonals */
  5105. -        A->me[k+2][k].re = A->me[k+2][k].im = 0.0;
  5106. -        }
  5107. -
  5108. -        /* test to see if matrix should split */
  5109. -        for ( k = k_min; k < k_max; k++ )
  5110. -        if ( zabs(A_me[k+1][k]) < MACHEPS*
  5111. -            (zabs(A_me[k][k])+zabs(A_me[k+1][k+1])) )
  5112. -        {
  5113. -            A_me[k+1][k].re = A_me[k+1][k].im = 0.0;
  5114. -            split = TRUE;
  5115. -        }
  5116. -
  5117. -    }
  5118. -    }
  5119. -    
  5120. -    /* polish up A by zeroing strictly lower triangular elements
  5121. -       and small sub-diagonal elements */
  5122. -    for ( i = 0; i < A->m; i++ )
  5123. -    for ( j = 0; j < i-1; j++ )
  5124. -        A_me[i][j].re = A_me[i][j].im = 0.0;
  5125. -    for ( i = 0; i < A->m - 1; i++ )
  5126. -    if ( zabs(A_me[i+1][i]) < MACHEPS*
  5127. -        (zabs(A_me[i][i])+zabs(A_me[i+1][i+1])) )
  5128. -        A_me[i+1][i].re = A_me[i+1][i].im = 0.0;
  5129. -
  5130. -    return A;
  5131. -}
  5132. -
  5133. -
  5134. -#if 0
  5135. -/* schur_vecs -- returns eigenvectors computed from the real Schur
  5136. -        decomposition of a matrix
  5137. -    -- T is the block upper triangular Schur matrix
  5138. -    -- Q is the orthognal matrix where A = Q.T.Q^T
  5139. -    -- if Q is null, the eigenvectors of T are returned
  5140. -    -- X_re is the real part of the matrix of eigenvectors,
  5141. -        and X_im is the imaginary part of the matrix.
  5142. -    -- X_re is returned */
  5143. -MAT    *schur_vecs(T,Q,X_re,X_im)
  5144. -MAT    *T, *Q, *X_re, *X_im;
  5145. -{
  5146. -    int    i, j, limit;
  5147. -    Real    t11_re, t11_im, t12, t21, t22_re, t22_im;
  5148. -    Real    l_re, l_im, det_re, det_im, invdet_re, invdet_im,
  5149. -        val1_re, val1_im, val2_re, val2_im,
  5150. -        tmp_val1_re, tmp_val1_im, tmp_val2_re, tmp_val2_im, **T_me;
  5151. -    Real    sum, diff, discrim, magdet, norm, scale;
  5152. -    static VEC    *tmp1_re=VNULL, *tmp1_im=VNULL,
  5153. -            *tmp2_re=VNULL, *tmp2_im=VNULL;
  5154. -
  5155. -    if ( ! T || ! X_re )
  5156. -        error(E_NULL,"schur_vecs");
  5157. -    if ( T->m != T->n || X_re->m != X_re->n ||
  5158. -        ( Q != MNULL && Q->m != Q->n ) ||
  5159. -        ( X_im != MNULL && X_im->m != X_im->n ) )
  5160. -        error(E_SQUARE,"schur_vecs");
  5161. -    if ( T->m != X_re->m ||
  5162. -        ( Q != MNULL && T->m != Q->m ) ||
  5163. -        ( X_im != MNULL && T->m != X_im->m ) )
  5164. -        error(E_SIZES,"schur_vecs");
  5165. -
  5166. -    tmp1_re = v_resize(tmp1_re,T->m);
  5167. -    tmp1_im = v_resize(tmp1_im,T->m);
  5168. -    tmp2_re = v_resize(tmp2_re,T->m);
  5169. -    tmp2_im = v_resize(tmp2_im,T->m);
  5170. -    MEM_STAT_REG(tmp1_re,TYPE_VEC);
  5171. -    MEM_STAT_REG(tmp1_im,TYPE_VEC);
  5172. -    MEM_STAT_REG(tmp2_re,TYPE_VEC);
  5173. -    MEM_STAT_REG(tmp2_im,TYPE_VEC);
  5174. -
  5175. -    T_me = T->me;
  5176. -    i = 0;
  5177. -    while ( i < T->m )
  5178. -    {
  5179. -        if ( i+1 < T->m && T->me[i+1][i] != 0.0 )
  5180. -        {    /* complex eigenvalue */
  5181. -        sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  5182. -        diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  5183. -        discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  5184. -        l_re = l_im = 0.0;
  5185. -        if ( discrim < 0.0 )
  5186. -        {    /* yes -- complex e-vals */
  5187. -            l_re = sum;
  5188. -            l_im = sqrt(-discrim);
  5189. -        }
  5190. -        else /* not correct Real Schur form */
  5191. -            error(E_RANGE,"schur_vecs");
  5192. -        }
  5193. -        else
  5194. -        {
  5195. -        l_re = T_me[i][i];
  5196. -        l_im = 0.0;
  5197. -        }
  5198. -
  5199. -        v_zero(tmp1_im);
  5200. -        v_rand(tmp1_re);
  5201. -        sv_mlt(MACHEPS,tmp1_re,tmp1_re);
  5202. -
  5203. -        /* solve (T-l.I)x = tmp1 */
  5204. -        limit = ( l_im != 0.0 ) ? i+1 : i;
  5205. -        /* printf("limit = %d\n",limit); */
  5206. -        for ( j = limit+1; j < T->m; j++ )
  5207. -        tmp1_re->ve[j] = 0.0;
  5208. -        j = limit;
  5209. -        while ( j >= 0 )
  5210. -        {
  5211. -        if ( j > 0 && T->me[j][j-1] != 0.0 )
  5212. -        {   /* 2 x 2 diagonal block */
  5213. -            /* printf("checkpoint A\n"); */
  5214. -            val1_re = tmp1_re->ve[j-1] -
  5215. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  5216. -            /* printf("checkpoint B\n"); */
  5217. -            val1_im = tmp1_im->ve[j-1] -
  5218. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  5219. -            /* printf("checkpoint C\n"); */
  5220. -            val2_re = tmp1_re->ve[j] -
  5221. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  5222. -            /* printf("checkpoint D\n"); */
  5223. -            val2_im = tmp1_im->ve[j] -
  5224. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  5225. -            /* printf("checkpoint E\n"); */
  5226. -            
  5227. -            t11_re = T_me[j-1][j-1] - l_re;
  5228. -            t11_im = - l_im;
  5229. -            t22_re = T_me[j][j] - l_re;
  5230. -            t22_im = - l_im;
  5231. -            t12 = T_me[j-1][j];
  5232. -            t21 = T_me[j][j-1];
  5233. -
  5234. -            scale =  fabs(T_me[j-1][j-1]) + fabs(T_me[j][j]) +
  5235. -            fabs(t12) + fabs(t21) + fabs(l_re) + fabs(l_im);
  5236. -
  5237. -            det_re = t11_re*t22_re - t11_im*t22_im - t12*t21;
  5238. -            det_im = t11_re*t22_im + t11_im*t22_re;
  5239. -            magdet = det_re*det_re+det_im*det_im;
  5240. -            if ( sqrt(magdet) < MACHEPS*scale )
  5241. -            {
  5242. -                det_re = MACHEPS*scale;
  5243. -            magdet = det_re*det_re+det_im*det_im;
  5244. -            }
  5245. -            invdet_re =   det_re/magdet;
  5246. -            invdet_im = - det_im/magdet;
  5247. -            tmp_val1_re = t22_re*val1_re-t22_im*val1_im-t12*val2_re;
  5248. -            tmp_val1_im = t22_im*val1_re+t22_re*val1_im-t12*val2_im;
  5249. -            tmp_val2_re = t11_re*val2_re-t11_im*val2_im-t21*val1_re;
  5250. -            tmp_val2_im = t11_im*val2_re+t11_re*val2_im-t21*val1_im;
  5251. -            tmp1_re->ve[j-1] = invdet_re*tmp_val1_re -
  5252. -                    invdet_im*tmp_val1_im;
  5253. -            tmp1_im->ve[j-1] = invdet_im*tmp_val1_re +
  5254. -                    invdet_re*tmp_val1_im;
  5255. -            tmp1_re->ve[j]   = invdet_re*tmp_val2_re -
  5256. -                    invdet_im*tmp_val2_im;
  5257. -            tmp1_im->ve[j]   = invdet_im*tmp_val2_re +
  5258. -                    invdet_re*tmp_val2_im;
  5259. -            j -= 2;
  5260. -            }
  5261. -            else
  5262. -        {
  5263. -            t11_re = T_me[j][j] - l_re;
  5264. -            t11_im = - l_im;
  5265. -            magdet = t11_re*t11_re + t11_im*t11_im;
  5266. -            scale = fabs(T_me[j][j]) + fabs(l_re);
  5267. -            if ( sqrt(magdet) < MACHEPS*scale )
  5268. -            {
  5269. -                t11_re = MACHEPS*scale;
  5270. -            magdet = t11_re*t11_re + t11_im*t11_im;
  5271. -            }
  5272. -            invdet_re =   t11_re/magdet;
  5273. -            invdet_im = - t11_im/magdet;
  5274. -            /* printf("checkpoint F\n"); */
  5275. -            val1_re = tmp1_re->ve[j] -
  5276. -              __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  5277. -            /* printf("checkpoint G\n"); */
  5278. -            val1_im = tmp1_im->ve[j] -
  5279. -              __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  5280. -            /* printf("checkpoint H\n"); */
  5281. -            tmp1_re->ve[j] = invdet_re*val1_re - invdet_im*val1_im;
  5282. -            tmp1_im->ve[j] = invdet_im*val1_re + invdet_re*val1_im;
  5283. -            j -= 1;
  5284. -        }
  5285. -        }
  5286. -
  5287. -        norm = v_norm_inf(tmp1_re) + v_norm_inf(tmp1_im);
  5288. -        sv_mlt(1/norm,tmp1_re,tmp1_re);
  5289. -        if ( l_im != 0.0 )
  5290. -        sv_mlt(1/norm,tmp1_im,tmp1_im);
  5291. -        mv_mlt(Q,tmp1_re,tmp2_re);
  5292. -        if ( l_im != 0.0 )
  5293. -        mv_mlt(Q,tmp1_im,tmp2_im);
  5294. -        if ( l_im != 0.0 )
  5295. -        norm = sqrt(in_prod(tmp2_re,tmp2_re)+in_prod(tmp2_im,tmp2_im));
  5296. -        else
  5297. -        norm = v_norm2(tmp2_re);
  5298. -        sv_mlt(1/norm,tmp2_re,tmp2_re);
  5299. -        if ( l_im != 0.0 )
  5300. -        sv_mlt(1/norm,tmp2_im,tmp2_im);
  5301. -
  5302. -        if ( l_im != 0.0 )
  5303. -        {
  5304. -        if ( ! X_im )
  5305. -        error(E_NULL,"schur_vecs");
  5306. -        set_col(X_re,i,tmp2_re);
  5307. -        set_col(X_im,i,tmp2_im);
  5308. -        sv_mlt(-1.0,tmp2_im,tmp2_im);
  5309. -        set_col(X_re,i+1,tmp2_re);
  5310. -        set_col(X_im,i+1,tmp2_im);
  5311. -        i += 2;
  5312. -        }
  5313. -        else
  5314. -        {
  5315. -        set_col(X_re,i,tmp2_re);
  5316. -        if ( X_im != MNULL )
  5317. -            set_col(X_im,i,tmp1_im);    /* zero vector */
  5318. -        i += 1;
  5319. -        }
  5320. -    }
  5321. -
  5322. -    return X_re;
  5323. -}
  5324. -
  5325. -#endif
  5326. //GO.SYSIN DD zschur.c
  5327.